Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

AIM: : To write a program for in MATLAB for 4th order approximation for second order derivative by following methods:

1. Central difference
2. Skewed right sided difference
3. Skewed left sided difference

OBJECTIVE:

• To evaluate the second order derivative of the function f(x) = exp(x)cos x, Using above mentioned order of approximations.
• To plot and compare the values and errors of the 3 methods derived.

EXPLANATION AND THEORY:

As the order of approximation get higher, we get higher accuracy. Hence, we are going to derive the three methods of approximation for an analytical function. The methods are following:

1.Central difference: Here the nodes selected are symmetrical i.e. the data is taken from both sides of point of focus. Using Taylor’s method, we get
∂2f/∂x2 = af(i − 2) + bf(i − 1) + cf(i) + df(i + 1) + ef(i + 2)
Here there are two nodes taken from both side of the point of focus i.e. f(i).
NOTE:- To select the number of nodes in consideration, we use following expression:
N=O + A -1 for central differencing
N = O + A for skewed differencing
Here N is no. of nodes, O is order of derivative and A is order of approximation.

Taylor’s table for Central difference:

f(X) dx f’(x) dx f’’(x) dx f’’’(x) dx f’’’’(x) dx f’’’’’(x)
af(i-1) a -2a 4a/2 -8a/6 16a/24 -32a/120
bf(i-2) b -b b/2 -b/6 b/24 -b/120
cf(i) c 0 0 0 0 0
df(i-1) d d d/2 d/6 d/24 d/120
ef(i-2) e 2e 4e/2 8e/6 16e/24 32e/120
af(i-2) +bf(i-1) +cf(i) +df(i+1) +ef(i+2) 0 0 1 0 0 ?

2.Skewed right sided difference: When the data point of the left nodes are not available (at initial boundary conditions) then we have to rely on the available data set.

∂2f/∂x2 = af(i) + bf(i + 1) + cf(i + 2) + df(i + 3) + ef(i + 4) + ff(i + 5)

Taylor’s table for Skewed right sided difference:

f(x) dx f’ (X) dx f’’ (x) dx f’’’(x) dx f’’’’(x) dx f’’’’’(x)

af(i) a 0 0 0 0 0
bf(i+1) b b b/2 b/6 b/24 b/120
cf(i+2) c 2c 4c/2 8c/6 16c/24 32c/120
df(i+3) d 3d 9d/2 27d/6 81d/24 243d/120
ef(i+4) e 4e 16e/2 64e/6 256e/24 1024e/120
ff(i+5) f 5f 25f/2 125f/6 625f/24 3125f/120
af(i) + bf(i+1) + cf(i+2) +df(i+3) + ef(i+4) +ff(i+5) 0 0 1 0 0 0

3.Skewed left sided difference: When the data point of the right nodes are not available (at boundary conditions) then we have to rely on the available data set

∂2f/∂x2 = af(i − 5) + bf(i − 4) + cf(i − 3) + df(i − 2) + ef(i − 1) + ff(i)

Taylor’s table for Skewed right sided difference:

f(x) dx f’ (X) dx f’’ (x) dx f’’’(x) dx f’’’’(x) dx f’’’’’(x)

af(i-5) a -5a 25a/2 -125a/6 625a/24 -3125a/120
bf(i-4) b -4b 16b/2 -64b/6 256b/24 -1024b/120
cf(i-3) c -3c 9c/2 -27c/6 81c/24 -243c/120
df(i-2) d -2d 4d/2 -8d/6 16d/24 -32d/120
ef(i-1) e -e e/2 -e/6 e/24 -e/120
ff(i) f 0 0 0 0 0
af(i-5) + bf(i-4) + cf(i-3) + df(i-2) + ef(i-1) +ff(i) 0 0 1 0 0 0

MATLAB CODE:

clear all
 close all
 clc

 %inputs

 x=pi/3;

 %stepsize range
 dx=linspace(pi/4,pi/400,2000);

 %Analytical derivative

 %f(x) = exp(x)*cos(x)
 %f'(x)=(exp(x)*cos(x)-(exp(x)*sin(x)=exp(x)*(cosx-sinx)
 %Double derivative f" = -2*exp(x)*sin(x)

 analytical_derivative = -2*exp(x)*sin(x);

 %Skewed right-sided(4th order approimation)
 %Solving for the coefficients
 A=[1 1 1 1 1 1;0 1 2 3 4 5; 0 1/2 4/2 9/2 16/2 25/2; 0 1/6 8/6 27/6 64/6 125/6; 0 1/24 16/24 81/24 256/24 625/24; 0 1/120 32/120 243/120 1024/120 3125/120];
 B=[0;0;1;0;0;0];
 x_r=inv(A)*B;

 %Skewed left-sided(4th order approimation)
 %Solving for the coefficients
 C=[1 1 1 1 1 1;0 -1 -2 -3 -4 -5; 0 1/2 4/2 9/2 16/2 25/2; 0 -1/6 -8/6 -27/6 -64/6 -125/6; 0 1/24 16/24 81/24 256/24 625/24; 0 -1/120 -32/120 -243/120 -1024/120 -3125/120];
 D=[0;0;1;0;0;0];
 x_l=inv(C)*D;

 %central differencing(4th order approimation)
 %Solving for the coefficients
 E=[1, 1, 1, 1, 1; -2, -1, 0, 1, 2; 4/2, 1/2, 0, 1/2, 4/2; -8/6, -1/6, 0, 1/6, 8/6; 16/24, 1/24, 0, 1/24, 16/24];
 F=[0;0;1;0;0];
 x_c=E^-1*F;

 %Anonymous Function

 f=@(x) exp(x).*cos(x);


 %numerical derivative

 right_difference=(x_r(1)*f(x)+x_r(2)*f(x+dx)+x_r(3)*f(x+2.*dx)+x_r(4)*f(x+3.*dx)+x_r(5)*f(x+4.*dx)+x_r(6)*f(x+5.*dx))./(dx.^2)

 left_difference=(x_r(1)*f(x)+x_r(2)*f(x-dx)+x_r(3)*f(x-2.*dx)+x_r(4)*f(x-3.*dx)+x_r(5)*f(x-4.*dx)+x_r(6)*f(x-5.*dx))./(dx.^2)

 cent_difference=(x_c(1)*f(x-2.*dx)+x_c(2)*f(x-1.*dx)+x_c(3)*f(x)+x_c(4)*f(x+dx)+x_c(5)*f(x+2.*dx))./(dx.^2)



 %absolute error

 error_r=abs(analytical_derivative-right_difference);
 error_l=abs(analytical_derivative-left_difference);
 error_c=abs(analytical_derivative-cent_difference);
 error=[error_r;error_l;error_c];



 figure(1)
 loglog(dx,error_r)
 grid on
 hold on
 loglog(dx,error_l)
 loglog(dx,error_c)
 ylabel('Absolute Error (Log Scale)')
 xlabel('Grid Spacing (dx in Log Scale)')
 legend('Right Differencing','Left Differencing','Central Differencing')

RESULTS AND OUTPUT:

CONCLUSION:

As we can see from the above figure the error for the central difference method is very less compared to the skewed methods.

Skewed right and left difference methods almost gives the same results in terms of error for a small value of time step.

• When data from the grid points is only available from either the left or right end i.e. at nodes closer to the boundaries,for computations that requires a high order approximations we make use of skewed schemes.Having a symmetric stencil for computation is not possioble in this scenarios like we do for central differencing schemes or having a low order accuracy.For better consistency of scheme to be maintained the values of dx are kept as low as possible for providing proper computational plots.
• A skewed scheme is useful when we have to predict using the boundary condition as we have only one direction to go. Hence the skewed scheme will be useful and we cannot use the central differencing scheme for such an application because it required information from both directions of the point.