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

Aim- Deriving a 4th-order approximation of a 2nd-order derivative using the Taylor Table method.

1. Right skewed difference scheme
2. Central difference scheme
3. Left Skewed difference scheme

Objective-

To evaluate the second order derivative of the function f(x)= exp(x)cos(x)

Using the above three order of approximations.

To plot and compare the values and errors of the 3 methods derived.

Theory and derivations: -

During solving a CFD problem we may come across a partial derivative equation which consists of 2nd order derivative. Since due to the complexity of the equation we may have to chose a numerical approach rather than analytical one. But when dealing with numerical approach we try to make the error margin as low as possible.

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 the following.

1. Right-Skewed Difference Scheme.

`(d^2f)/dx^2 =af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)`

Taylor table for right-sided skewed difference scheme

Stensils	f(i)	dxf(i)	dxf''(i)	dxf'''(i)	dxf''''(i)	dx'''''f(i)
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
gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120
af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)	0	0	1	0	0	0

table transform into the matrix of type AX=B

`[[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]]*[(a),(b),(c),(d),(e),(f)] = [(0),(0),(1),(0),(0),(0)]`

For further solving the values of unknown

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.

2. Central Difference Scheme.

`(d^2f)/dx^2 =af(i-1)+bf(i-2)+cf(i)+df(i+1)+ef(i+2)`

Taylor table for central difference scheme

Stensils	f(i)	dxf'(i)	dxf''(i)	dxf'''(i)	dxf''''(i)	dxf'''''(i)
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-1)+bf(i-2)+cf(i)+df(i+1)+ef(i+2)	0	0	1	0	0	?

table transform into the matrix of type AX=B

[[1,1,1,1,1],[-2,-1,0,1,2],[4/2,1/2,0,1/2,4/2],[1,2,3,4,5],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]]*[(a),(b),(c),(d),(e),(f)] = [(0),(0),(1),(0),(0),(0)]

For further solving the values of unknown

Values = a = -0.0833 b= -1.3333 c= -2.5000

d= 1.3333 e= -0.0833

3. Left-Skewed Difference Scheme.

`(d^2f)/dx^2 =af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)`

Taylor table for Left skewed difference scheme

Stensils	f(i)	dxf'(i)	dxf''(i)	dxf'''(i)	dxf''''(i)	dxf'''''(i)
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
gf(i)	g	0	0	0	0	0
af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)	0	0	1	0	0	?

table transform into the matrix of type AX=B

[[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]]*[(a),(b),(c),(d),(e),(f)] = [(0),(0),(1),(0),(0),(0)]

For further solving the values of unknown

Values = a = -0.0833 b= 5.0833 c= 13.0000

d= 17.8333 e= -12.8333 f= 3.7500

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 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 '\')

cf1=Q\D;

%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 '\')

cf2=P\C;

%for left 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 right skewed difference scheme
    right_sdiff(i)=(cf1(1)*f(x)+cf1(2)*f(x+dx(i))+cf1(3)*f(x+2*dx(i))+cf1(4)*f(x+3*dx(i))+cf1(5)*f(x+4*dx(i))+cf1(6)*f(x+5*dx(i)))/(dx(i))^2;
    err_right(i)=abs(right_sdiff(i)-d2f);
    
    %calculating value and error by central difference scheme
    central_diff(i)=(cf2(1)*f(x-2*dx(i))+cf2(2)*f(x-dx(i))+cf2(3)*f(x)+cf2(4)*f(x+dx(i))+cf2(5)*f(x+2*dx(i)))/(dx(i))^2;
    err_central(i)=abs(central_diff(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 ');

Conclusion-

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