Formulation of 4th order central and skewed approximation using taylor table and implementing in matlab

Taylor series is used to aproximate the values of derivatives at a certain through various approximate scheme. The drivative at a node is measured with the help of the data at the surrounding node.The order of the scheme is indicates that the how fast the error decreases when the node spacing is reducing. The higher the order, the more accurate the scheme, but also it increases the computational cost. Derivation of central, left skewed and right skewed are given below for the fourth order approximation of second order derivative.

1. Fourth order central difference of the second order derivative:

`frac{partial^2u}{partial x^2}= frac{au(i-2)+bu(i-1)+cu(i)+du(i+1)+eu(i+2)}{Deltax^2}`

in central difference we use the A+O-1=4+2-1=5 to calculate the no of terms.

in skewed schemes, we represent the A+O=6 to calculate the no of terms.

where A represent the accuracy and O represent the order.

Formulation of Taylor table:

1st term=f(i)

2nd term=`fprime Deltax`

3rd term= `fprimeprimeDeltax^2`

4th term=`fprimeprimeprimeDeltax^3`

5th term= `fprimeprimeprimeprime Deltax^4`

6th term= `fprimeprimeprimeprimeprime Deltax^5`

Tatlor series Stencils 1st term 2nd term 3rd term 4th term 5th term

au(i-2) a -2a 2a -4a/3 2a/3

bu(i-1) b -b b/2 -b/6 b/24

cu(i) c 0 0 0 0

du(i+1) d d d/2 d/6 d/24

eu(i+2) e 2e 2e 4e/3 2e/3

sum 0 0 1 0 0

Basically taylor table is made by expanding the well known taylor series.The values of a,b,c,d e are solved by using the matlab code which is given below

clear all
close all
clc
%for central difference
    A=[ 1 1 1 1 1;
    -2 -1 0 1 2;
    2 1/2 0 1/2 2;
    -4/3 -1/6 0 1/6 4/3;
    2/3 1/24 0 1/24 2/3]
B= [0; 0; 1; 0; 0]
solution= inv(A)*B
check=sum(solution)

And output obtained from this code is

a=-0.0833
b=1.3333
c=-2.5000
d=1.3333
e=-0.0833

2. Fourth order skewed right sided difference of the second order derivative

`frac{partial^2u}{partialx}=frac{au(i)+bu(i+1)+cu(i+2)+du(i+3)+eu(i+4)+fu(i+5)}{Deltax^2}`

Formulation of Taylor table

Taylor series stencils 1st term 2nd term 3rd term 4th term 5th term 6th term

au(i) a 0 0 0 0 0

bu(i+1) b b b /2 b/6 b/24 b/120

cu(i+2) c 2c 2c 4c/3 2c/3 32c/120

du(i+3) d 3d 9d/2 27d/6 81d/24 243d/120

eu(i+4) e 4e 8e 64e/6 256e/24 1024e/120

fu(i+5) f 5f 25f/2 125f/6 625f/24 3125f/120

sum 0 0 1 0 0

values of a,b,c,d,e,f are solved by using the matlab code which is given below

clear all
close all
clc
%for right skewed difference
A=[ 1 1 1 1 1 1;
    0 1 2 3 4 5;
    0 1/2 2 9/2 8 25/2;
    0 1/6 4/3 27/6 64/6 125/6;
    0 1/24 2/3 81/24 256/24 625/24;
    0 1/120 32/120 243/120 1024/120 3125/120]
B= [0; 0; 1; 0; 0;0]
solution= inv(A)*B
check=sum(solution)

And the output obtained from above code is

a=3.7500
b=-12.8333
c=17.8333
d=-13.0000
e=5.0833
f=-0.8333

3.Fourth order skewed left sided difference of the second order derivative

`frac{partial^2u}{partialx}=frac{au(i-5)+bu(i-4)+cu(i-3)+du(i-2)+eu(i-1)+fu(i)}{Deltax^2}`

Formulation of Taylor table

Taylor series stencils 1st term 2nd term 3rd term 4th term 5th term 6th term

au(i) a 0 0 0 0 0

bu(i-1) b -b b/2 -b/6 b/24 -b/120

cu(i-2) c -2c 2c -4c/3 2c/3 -32c/120

du(i-3) d -3d 9d/2 -27d/6 81d/24 -243d/120

eu(i-4) e -4e 8e -64e/6 256e/24 -1024e/120

fu(i-5) f -5f 25f/2 -125f/6 625f/24 -3125f/120

sum 0 0 1 0 0 0

values of a,b,c,d,e,f are solved by using matlab code which is given below

clear all
close all
clc
%for left skewed difference
A=[ 1 1 1 1 1 1;
    0 -1 -2 -3 -4 -5;
    0 1/2 2 9/2 8 25/2;
    0 (-1/6) (-4/3) (-27/6) (-64/6) (-125/6);
    0 1/24 2/3 81/24 256/24 625/24;
    0 (-1/120) (-32/120) (-243/120) (-1024/120) (-3125/120)]
B= [0; 0; 1; 0; 0;0]
solution= inv(A)*B
check=sum(solution)

And the output obtained from the above code is

a=3.7500
b=-12.8333
c=17.8333
d=-13.0000
e=5.0833
f=-0.8333

Matlab code

%fourth order approximation
%given function
%f(x)= (exp(x)*cos(x));
x=pi/3;
dx=linspace(pi/40,pi/400,20);
analytical_derivative= (-2)*exp(x)*sin(x);
for i=1:length(dx)
    central_differencing(i)=((-(1/12)*epxcs(x-(2*dx(i))))+((4/3)*epxcs(x-dx(i)))-(2.5000*epxcs(x))+((4/3)*epxcs(x+dx(i)))-((1/12)*epxcs(x+(2*dx(i)))))/(dx(i)^2);
    skewed_left_differencing(i)=(((15/4)*epxcs(x))-((77/6)*epxcs(x-dx(i)))+((107/6)*epxcs(x-(2*dx(i))))-(13*epxcs(x-(3*dx(i))))+((61/12)*epxcs(x-(4*dx(i))))-((5/6)*epxcs(x-(5*dx(i)))))/(dx(i)^2);
    skewed_right_differencing(i)=(((15/4)*epxcs(x))-((77/6)*epxcs(x+dx(i)))+((107/6)*epxcs(x+(2*dx(i))))-(13*epxcs(x+(3*dx(i))))+((61/12)*epxcs(x+(4*dx(i))))-((5/6)*epxcs(x+(5*dx(i)))))/(dx(i)^2);
end
%error calculation
central_error=abs(analytical_derivative-central_differencing);
skewed_left_error=abs(analytical_derivative-skewed_left_differencing);
skewed_right_error=abs(analytical_derivative-skewed_right_differencing);
loglog(dx, central_error, \'--r\', \'linewidth\',2)
grid on
hold on
loglog(dx, skewed_left_error, \'--b\', \'linewidth\',2)
hold on
loglog(dx, skewed_right_error, \'--g\', \'linewidth\',2)
legend(\'central error\', \'skewed_left_error\', \'skewed_right_error\')
title(\'comparison of error from different scheme\')
xlabel(\'dx values\')
ylabel(\'Error values\')

Function code which is used in above program

function out= epxcs(x)
out= exp(x)*cos(x);
end

Results discussion:

From the above graph, as we observe that the error of central difference is lower as compared to the skewed schemes. This is because the quality of approximation improves as the stencils point get close to the node x. In central difference scheme the farthest point of the stencil is at a distance of 2h , where for the skewed stencil, the farthest point are at a distance of 5h. Also the consistency of central difference scheme is much batter than the skewed schemes.

However, skewed schemes are useful at a node closer to the boundaries of the domain, where the symmetric stencil is either not possible. in skewed schemes we get the information of nodes at one side either left or right that\'s why error grow, while in central difference schemes we get the informations of nodes on both sides that\'s why error reduce. suppose we are at the start of grid, there we donot have a information from the previous grid, in these types of problems skewed schemes are helpful.