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

Aim: - To derive and write a program for in MATLAB for 4^th order approximation for second order derivative by following methods-

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

Objective: -

1. To estimate the alternate order outgrowth of the function
2.To compass and compare the values and crimes of the 3 styles deduced.

Theory and derivations: -

During working a CFD problem we may come across partial secondary equation which consists 2nd order outgrowth. Since due to complexity of the equation we may have to chose numerical approach rather than logical. But when dealing with numerical approach we try to make the error periphery as low as possible.

As the order of approximation get advanced, we get advanced delicacy. Hence, we're going to decide the three styles of approximation for an logical function. The styles are following.

1.Central difference- Then the bumps named are symmetrical i.e. the data is taken from both sides of point of focus. Using Taylor’s system, we get

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=p + q -1 for central differencing

N = p + q for skewed differencing

Here N is no. of nodes, p is order of derivative and q is order of approximation.)

Taylor’s table for Central difference:

	f(X)	dx f’(x)	dx2 f’’(x)	dx3 f’’’(x)	dx4 f’’’’(x)	dx5 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	?

This table transform into matrix of type AX=B

i.e. Values = a = -0.0833 b= 1.3333 c= -2.5000 d= 1.3333 e= -0.0833

These coefficients will further used in code to get the value of derivative.

1. 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.

Taylor’s table for Skewed right sided difference:

	f(x)	dx f’ (X)	dx2 f’’ (x)	dx3 f’’’(x)	dx4f’’’’(x)	dx5f’’’’’(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

This table transform into matrix of type AX=B:[1111110123450122921622520168627664612560124162481242562462524011203212024312010241203125120]">

unknown i.e. Values = a = 3.7500 b= -12.8333 c= 17.8333

d= -13.0000 e= 5.0833 f= -0.8333

These coefficients will further used in code to get the value of derivative.

1. 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.

Taylor’s table for Skewed right sided difference:

	f(x)	dx f’ (X)	dx2 f’’ (x)	dx3 f’’’(x)	dx4f’’’’(x)	dx5f’’’’’(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

This table transform into matrix of type AX=B:

unknown i.e. Values = a = -0.8333 b= 5.0833 c=13.0000

d= 17.8333 e= -12.8333 f= 3.7500

These coefficients will further used in code to get the value of derivative.

MATLAB code: -

clear all
close all
clc

%analytical function defination and value assignment
x=pi/3;
dx=linspace(pi/10,pi/1000,100);
f=@(x)(exp(x)*cos(x));
g=@(x)(-2*exp(x)*sin(x));
d2f=g(x);
func=f(x);



%for central difference scheme 
%coefficients taken from taylor's table
P = [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];
    
   
    
    
C=[0;0;1;0;0];
%as P*a=C i.e. a=inverse(P)*c(fastest way to getby '\')

cf1=P\C;

%for right skewwd differnce scheme
Q= [1 1 1 1 1 1
    0 1 2 3 4 5
    0 1/2 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];
%as Q*b=s i.e. b=inverse(Q)*s(fastest way to getby '\')

cf2=Q\D;

%for right skewwd differnce scheme
R= [1 1 1 1 1 1
    -5 -4 -3 -2 -1 0
    25/2 16/2 9/2 2 1/2 0
    -125/6 -64/6 -27/6 -8/6 -1/6 0
    625/24 256/24 81/24 16/24 1/24 0
    -3125/120 -1024/120 -243/120 -32/120 -1/120 0];
E=[0;0;1;0;0;0];
%as n*b=s i.e. b=inverse(n)*s(fastest way to getby '\')

cf3=R\E;

%calculating values of derivative using CDS RSDC,LSDC
%defining variable size befor using in loop
central_diff=zeros(1,100);
err_central=zeros(1,100);
right_sdiff=zeros(1,100);
err_right=zeros(1,100);
left_sdiff=zeros(1,100);
err_left=zeros(1,100);

for i=1:length(dx)
    %calculating value and error by central difference scheme
    central_diff(i)=(cf1(1)*f(x-2*dx(i))+cf1(2)*f(x-dx(i))+cf1(3)*f(x)+cf1(4)*f(x+dx(i))+cf1(5)*f(x+2*dx(i)))/(dx(i))^2;
    err_central(i)=abs(central_diff(i)-d2f);
    
    %calculating value and error by right skewed difference scheme
    right_sdiff(i)=(cf2(1)*f(x)+cf2(2)*f(x+dx(i))+cf2(3)*f(x+2*dx(i))+cf2(4)*f(x+3*dx(i))+cf2(5)*f(x+4*dx(i))+cf2(6)*f(x+5*dx(i)))/(dx(i))^2;
    err_right(i)=abs(right_sdiff(i)-d2f);
    
    %calculating value and error by left skewed difference scheme
   left_sdiff(i)=(cf3(1)*f(x-5*dx(i))+cf3(2)*f(x-4*dx(i))+cf3(3)*f(x-3*dx(i))+cf3(4)*f(x-2*dx(i))+cf3(5)*f(x-dx(i))+cf3(6)*f(x))/(dx(i))^2;
   err_left(i)=abs(left_sdiff(i)-d2f); 

end

%plotting values of derivative
subplot(2,1,1)
hold on
fplot(d2f,'--')
xlim([-0.05 0.35]);
plot(dx,central_diff,'linewidth',2,'color','r');
plot(dx,right_sdiff,'linewidth',2,'color','b');
plot(dx,left_sdiff,'linewidth',2,'color','g');
legend('analytical value at x=pi/3','central diff.','right skewed diff.','left skewed diff.','location','best')
xlabel('Grid siza(dx)');
ylabel('Value');
title('Comparision b\w value of derivative calculated by diff. order of approximation') 
hold off

%plotting errors in approximation of derivative
subplot(2,1,2)
loglog(dx,err_central,'linewidth',2,'color','r');
hold on
xlim([.002 0.4]);
loglog(dx,err_right,'linewidth',2,'color','b');
loglog(dx,err_left,'linewidth',2,'color','g');
legend('central diff.','right skewed diff.','left skewed diff.','location','best')
xlabel('Grid siza(dx)');
ylabel('Absolute error');
title('Comparision b/w 4th order approximation method for error in calculating d2f/dx2 ');

steps involved: -

1. We assign the value of x and range of dx
2. We define an anonymous function of logical function and its outgrowth.
3. Use the matrices from above bandied section and answered it using inversion system i.e.
   layoff = B == > X = A-1 B
4. And stored the portions in cf1, cf2 and cf3.
5. Assigned the rows and column to the applicable matrices.
6. Loop the order of approximation expression in a for circle under the range of dx. Then absolute error is also calculated from the secondary value calculated analytically.
7. The calculated valued are colluded to do a better comparision.

Result:-

1. In first plot we can see that the value of all the approximation approached toward the analytically calculated value as the size of dx gets smaller.
2. On the second plot we can see that the error in the central difference is lowest with respect to other two.