NUMERICAL DIFFERENTIATION AND INTEGRATION

 

OBJECTIVE

  • To perform  numerical differentiation and integration

LEARNING OUTCOMES

  • After the completion of this experiment students will be able to Implement numerical integration and differentiation.

SOFTWARE USED:

                 MATLAB® R2017a

THEORY

The first derivative of the function f(x), which we write as df /dx , is the slope of the tangent line to the function at the point x. To put this in non-graphical terms, the first derivative tells us how whether a function is increasing or decreasing, and by how much it is increasing or decreasing .

The second derivative of a function is the derivative of the derivative of that function. We write it as  d2 f /dx2 . While the first derivative can tell us if the function is increasing or decreasing, the second derivative tells us if the first derivative is increasing or decreasing. If the second derivative is positive, then the first derivative is increasing, so that the slope of the tangent line to the function is increasing as x increases

sin – Sine of argument in radians

cos - cos of argument in radians




Simpsons rule

Formula:  (h/3)*[(y0+yn)+2*(y3+y5+..odd term)+4*(y2+y4+y6+...even terms)]

h= (b-a)/n

b is upper limit, a is lower limit and n is number of sub intervals. y0 and yn are first and last term.

 

1. Realize the functions sin t, cos t, sinht and cosht for the vector t = [-0, 10] with increment 0:01

PROGRAM

t=-0:0.01:10;

a=sin (t);

b=cos (t);

c= sinh (t);

d=cosh (t);

subplot(2,2,1)

plot(t,a)

subplot(2,2,2)

plot(t,b)

subplot(2,2,3)

plot(t,c)

subplot(2,2,4)

plot(t,d)

 OUTPUT


2.   Compute the first and second derivatives of these functions using built in tools such as grad 

PROGRAM

t=0:0.01:10;

a=sin (t);

subplot(3,3,1)

plot(t,a)

title('sint')

b=cos (t);

c= sinh (t);

d=cosh (t);

da=gradient(a);

subplot(3,3,2)

plot(t,da)

title('first derivative of sint')

d2a=gradient(da);

subplot(3,3,3)

plot(t,d2a)

title('second derivative of sint')

db=gradient(b);

subplot(3,3,4)

plot(t,db)

title('first derivative of cost')

d2b=gradient(db);

subplot(3,3,5)

plot(t,d2b)

title('second  derivative of cost')

dc=gradient(c);

subplot(3,3,6)

plot(t,dc)

title('first derivative of sinht')

d2c=gradient(dc);

subplot(3,3,7)

plot(t,d2c)

title('second derivative of sinht')

dd=gradient(d);

subplot(3,3,8)

plot(t,dd)

title('first derivative of cosht')

d2d=gradient(dd);

subplot(3,3,9)

plot(t,d2d)

title('second  derivative of cosht')

OUTPUT


3. Compute the first and second derivatives of sint, cost, sinht, cosht functions using built in tools such as grad and plot the derivatives over the respective functions for the vector

 t = [-5, 5] with increment 0:01

 PROGRAM

t=-5:0.01:5;

a=sin (t);

subplot(3,3,1)

plot(t,a)

title('sint')

b=cos (t);

c= sinh (t);

d=cosh (t);

da=gradient(a);

subplot(3,3,2)

plot(t,da)

title('first derivative of sint')

d2a=gradient(da);

subplot(3,3,3)

plot(t,d2a)

title('second derivative of sint')

db=gradient(b);

subplot(3,3,4)

plot(t,db)

title('first derivative of cost')

d2b=gradient(db);

subplot(3,3,5)

plot(t,d2b)

title('second  derivative of cost')

dc=gradient(c);

subplot(3,3,6)

plot(t,dc)

title('first derivative of sinht')

d2c=gradient(dc);

subplot(3,3,7)

plot(t,d2c)

title('second derivative of sinht')

dd=gradient(d);

subplot(3,3,8)

plot(t,dd)

title('first derivative of cosht')

d2d=gradient(dd);

subplot(3,3,9)

plot(t,d2d)

title('second  derivative of cosht')

 OUTPUT


4. Familiarise numerical integration tools used in matlab

a. Create the function f(x)=ex2(lnx)2. Evaluate integral  from 0 to infinity

PROGRAM

fun = @(x) exp(-x.^2).*log(x).^2;

q = integral(fun,0,Inf)

OUTPUT

q = 1.9475

b. Create the function f(x)=1/(x3−2xc) with one parameter, c. Evaluate the integral from x=0 to x=2 at c=5.

PROGRAM

fun = @(x,c) 1./(x.^3-2*x-c);

q = integral(@(x) fun(x,5),0,2)

OUTPUT

q = -0.4605

c. Create the function f(x)=ln(x). Evaluate the integral from x=0 to x=1

PROGRAM

               fun = @(x)log(x);

               q1 = integral(fun,0,1)
OUTPUT
               q1 = -1.000000
 

5. Realize the function f(t) =4t2 +3 and plot it for the vector [-5,5] with increment 0.01

t=-5:0.01:5;

y=4*(t.^2)+3;

plot(t,y,'k','linewidth',2)

xlabel 't'

ylabel 'f(t)'

OUTPUT

6. Use general integration tool to compute 



General integration tool

PROGRAM

clc;

clear all;

fun = @(t) t + 2;

q = integral(fun,-2,2)

OUTPUT

 q=8

7. Repeat using trapezoidal method and Simpsons method

 PROGRAM

t=-2:.5:2;

y=t+2;

q=trapz(t,y)

 OUTPUT

 q=8


Trapezoidal method using loop

PROGRAM

clc;

clear all;

f=@(x)x+2; %Change here for different function

 a=input('Enter lower limit a: '); % exmple a=1

 b=input('Enter upper limit b: ');  % exmple b=2

 n=input('Enter the no. of subinterval: ');  % exmple n=10

 h=(b-a)/n;

 sum=0;

for k=1:1:n-1

  x(k)=a+k*h;

  y(k)=f(x(k));

  sum=sum+y(k);

end

% Formula:  (h/2)*[(y0+yn)+2*(y2+y3+..+yn-1)]

answer=h/2*(f(a)+f(b)+2*sum);

fprintf('\n The value of integration is %f',answer);

OUTPUT 

answer=8

Simpson method using loop

PROGRAM

clc;

clear all;

f=@(x)x+2; %Change here for different function

 a=input('Enter lower limit a: '); % exmple a=-2

 b=input('Enter upper limit b: ');  % exmple b=2

 n=input('Enter the  number of sub-intervals n: '); % exmple n=16

 h=(b-a)/n;

 

for k=1:1:n

  x(k)=a+k*h;

  y(k)=f(x(k));

end

 so=0;se=0;

for k=1:1:n-1

    if rem(k,2)==1

       so=so+y(k);%sum of odd terms

     else

       se=se+y(k); %sum of even terms

    end

end

% Formula:  (h/3)*[(y0+yn)+4*(y3+y5+..odd term)+2*(y2+y4+y6+...even terms)]

answer=h/3*(f(a)+f(b)+4*so+2*se);

fprintf('\n The value of integration is %f',answer);

 

8. Compute     


 
 using above three methods

PROGRAM

i)  t=0:.5:1000;

yi=@(t)exp((-t.^2)/2);

qi = integral(yi,0,1000)

answeri=(1/sqrt(2*pi))*qi

OUTPUT

answeri= 0.500

PROGRAM

ii)   t=0:.5:1000;

y=exp((-t.^2)/2);

q=trapz(t,y);

answer= (1/sqrt(2*pi))*q

OUTPUT

answer= 0.500

PROGRAM

iii)

 clc;

clear all;

f=@(x)exp((-x.^2)/2); %Change here for different function

 a=input('Enter lower limit a: '); % exmple a=1

 b=input('Enter upper limit b: ');  % exmple b=2

 n=input('Enter the  number of sub-intervals n: '); % exmple n=16

 h=(b-a)/n;

 

for k=1:1:n

  x(k)=a+k*h;

  y(k)=f(x(k));

end

 so=0;se=0;

for k=1:1:n-1

              if rem(k,2)==1

                 so=so+y(k);%sum of odd terms

              else

               se=se+y(k); %sum of even terms

             end

end

% Formula:  (h/3)*[(y0+yn)+2*(y3+y5+..odd term)+4*(y2+y4+y6+...even terms)]

q=h/3*(f(a)+f(b)+4*so+2*se);

answer= (1/sqrt(2*pi))*q

fprintf('\n The value of integration is %f',answer);

 

OUTPUT

Enter lower limit a: 0

Enter upper limit b: 1000

Enter the  number of sub-intervals n: 1000

answer =0.4976

 The value of integration is 0.497603

Comments

Popular posts from this blog

Coin Toss and Level crossing problem