Fourth Order Approximation of Second Derivative of a function using Central Right sided and Left sided Difference Shemes

Fourth Order Approximation of Second Derivative of a Function

1. Central Difference Scheme

In Central difference scheme the number of nodes required to get fourth order approximation for a second derivative of a function is =5 ( A+O-1) where A is order of approximation desired and O is the order of the derivative. In this case A=4 and O=2.

Let us assume our hypothesis equation to be:

`(/_\\ x)^2 f\'\'(i)=af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`...(1)

We will have to expand terms of the RHS of the equation (1) that is f(i-2),f(i-1),f(i+1),f(i+2) about f(i) using taylor series.

`f(i-2)=f(i)-2/_\\xf\'(i)+((2/_\\x)^2f\'\'(i))/(2!)-((2/_\\x)^3f\'\'\'(i))/(3!)+((2/_\\x)^4f\'\'\'\'(i))/(4!)...`

`f(i-1)=f(i)-/_\\xf\'(i)+((/_\\x)^2f\'\'(i))/(2!)-((/_\\x)^3f\'\'\'(i))/(3!)+((/_\\x)^4f\'\'\'\'(i))/(4!)...`

`f(i+1)=f(i)+/_\\xf\'(i)+((/_\\x)^2f\'\'(i))/(2!)+((/_\\x)^3f\'\'\'(i))/(3!)+((/_\\x)^4f\'\'\'\'(i))/(4!)...`

`f(i+2)=f(i)+2/_\\xf\'(i)+((2/_\\x)^2f\'\'(i))/(2!)+((2/_\\x)^3f\'\'\'(i))/(3!)+((2/_\\x)^4f\'\'\'\'(i))/(4!)...`

Creating a Taylor\'s Table to evaluate the coefficients a,b,c,d and e for the equation (1) using the above Taylor\'s series expansions of f(i-2),f(i-1),f(i+1),f(i+2).

Each row of the Taylor\'s table is obtained from the coefficients of the Taylor\'s series on multiplying a with f(i-2) expansion,b with f(i-1) expansion,c with f(i),d with f(i+1) expansion and e with f(i+2) expansion.

In the Taylors table we have equated the sum of coefficients of each of the f(i), f\'(i), f\'\'\'(i), f\'\'\'\'(i) to zeros and f\'\'(i) to 1 because we desire to obtain an expression in terms of f\'\'(i) in the RHS of the hypothesis equation (1).

The last column of the table is not required because we need 5 equations to solve for 5 unknowns but it is shown just to show the process of filling the table.

The linear system of equations obtained from Taylor\'s table are:

a+b+c+d+e=0

-2a-b-d+2e=0

2a+b/2+d/2+2e=1

-8a/6-b/6+d/6+8e/6=0

16a/24+b/24+d/24+16e/24=0

We can solve the system of equations using Matlab and get the value of the coefficients:

%% Computing the value of Co-eff for Central difference Scheme

A_CD=[1,1,1,1,1;-2,-1,0,1,2;2,1/2,0,1/2,2;-8/6,-1/6,0,1/6,8/6;16/24,1/24,0,1/24,16/24];
B=[0;0;1;0;0];
X_CD=A_CD\\B;

The above matlab code computes for the coefficients values(a,b,c,d,e) in a matrix X_CD. The A_CD is the coefficients Matrix of the above system of linear equations, B matrix gives the residuals of the linear equations. The equation is of the form A_CD*X_CD=B.

So, X_CD= (Inversion of A_CD matrix) * B

The output from Matlab is:

So the hypothesis equation for CDS scheme can be approximately written as:

`(/_\\ x)^2 f\'\'(i)=-f(i-2)/12+(4f(i-1))/3-(5f(i))/2+(4f(i+1))/3-f(i+2)/12+O(/_\\x)^6`

`f\'\'(i)=(-f(i-2)/12+(4f(i-1))/3-(5f(i))/2+(4f(i+1))/3-f(i+2)/12)/(/_\\x)^2+O(/_\\x)^4`

`f\'\'(i)~~(-f(i-2)+16f(i-1)-30f(i)+16f(i+1)-f(i+2))/(12(/_\\x)^2)`...(2)

We will use this equation (2) in Matlab program to compute Central Difference Scheme fourth order approximation of second derivative and corresponding error from exact second derivative.

2. Skewed Right Sided Difference Scheme

In skewed right sided difference scheme the number of nodes required to get fourth order approximation for a second derivative of a function is =6 ( A+O) where A is order of approximation desired and O is the order of the derivative. In this case A=4 and O=2.

Let us assume our hypothesis equation to be:

`(/_\\ x)^2 f\'\'(i)=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)` ....(3)

We will have to expand terms of the RHS of the equation (3) that is f(i+1), f(i+2), f(i+3), f(i+4) and f(i+5) about f(i) using taylor series.

`f(i+1)=f(i)+/_\\xf\'(i)+((/_\\x)^2f\'\'(i))/(2!)+((/_\\x)^3f\'\'\'(i))/(3!)+((/_\\x)^4f\'\'\'\'(i))/(4!)+((/_\\x)^5f\'\'\'\'(i))/(5!)+...`

`f(i+2)=f(i)+2/_\\xf\'(i)+((2/_\\x)^2f\'\'(i))/(2!)+((2/_\\x)^3f\'\'\'(i))/(3!)+((2/_\\x)^4f\'\'\'\'(i))/(4!)+((2/_\\x)^5f\'\'\'\'\'(i))/(5!)+...`

`f(i+3)=f(i)+3/_\\xf\'(i)+((3/_\\x)^2f\'\'(i))/(2!)+((3/_\\x)^3f\'\'\'(i))/(3!)+((3/_\\x)^4f\'\'\'\'(i))/(4!)+((3/_\\x)^5f\'\'\'\'\'(i))/(5!)...`

`f(i+4)=f(i)+4/_\\xf\'(i)+((4/_\\x)^2f\'\'(i))/(2!)+((4/_\\x)^3f\'\'\'(i))/(3!)+((4/_\\x)^4f\'\'\'\'(i))/(4!)+((4/_\\x)^5f\'\'\'\'\'(i))/(5!)...`

`f(i+5)=f(i)+5/_\\xf\'(i)+((5/_\\x)^2f\'\'(i))/(2!)+((5/_\\x)^3f\'\'\'(i))/(3!)+((5/_\\x)^4f\'\'\'\'(i))/(4!)+((5/_\\x)^5f\'\'\'\'\'(i))/(5!)...`

Creating a Taylor\'s Table to evaluate the coefficients a,b,c,d,e and g for the equation (3) using the above Taylor\'s series expansions of f(i+1),f(i+2),f(i+3),f(i+4),f(i+5).

Each row of the Taylor\'s table is obtained from the coefficients of the Taylor\'s series on multiplying a with f(i), b with f(i+1) expansion, c with f(i+2) expansion, d with f(i+3) expansion, e with f(i+4) expansion and g with f(i+5) expansion.

In the Taylors table we have equated the sum of coefficients of each of the f(i), f\'(i), f\'\'\'(i), f\'\'\'\'(i), f\'\'\'\'\'(i) to zeros and f\'\'(i) to 1 because we desire to obtain an expression in terms of f\'\'(i) in the RHS of the hypothesis equation (3).

The last column of the table is not required because we need 6 equations to solve for 6 unknowns but it is shown just to show the process of filling the table.

The linear system of equations obtained from Taylor\'s table are:

a+b+c+d+e+g=0

b+2c+3d+4e+5g=0

b/2+2c+9/2d+8e+25g/2=1

b/6+8c/6+27d/6+64e/6+125g/6=0

b/24+16c/24+81d/24+256e/24+625g/24=0

b/120+32c/120+243d/120+1024e/120+3125g/120=0

We can solve the system of equations using Matlab and get the value of the coefficients:

%% Computing the value of Co-eff for Skewed Right Sided Scheme
B=[0;0;1;0;0;0];
A_R=[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];
X_R=A_R\\B;

The above matlab code computes for the coefficients values(a,b,c,d,e,g) in a matrix X_R. The A_R is the coefficients Matrix of the above system of linear equations, B matrix gives the residuals of the linear equations. The equation is of the form A_R*X_R=B.

So, X_R= (Inversion of A_R matrix) * B

The output from Matlab is:

So the hypothesis equation for RDS scheme can be approximately written as:

`(/_\\ x)^2 f\'\'(i)=(15f(i))/4-(77f(i+1))/6+(107f(i+2))/6-13f(i+3)+(61f(i+4))/12-(5f(i+5))/6+O(/_\\x)^6`

`f\'\'(i)=((15f(i))/4-(77f(i+1))/6+(107f(i+2))/6-13f(i+3)+(61f(i+4))/12-(5f(i+5))/6)/(/_\\x)^2+O(/_\\x)^4`

`f\'\'(i)~~(45f(i)-154f(i+1)+214f(i+2)-156f(i+3)+61f(i+4)-10f(i+5))/(12(/_\\x)^2)`...(4)

We will use this equation (4) in Matlab program to compute Skewed Right sided Difference Scheme fourth order approximation of second derivative and corresponding error from exact second derivative.

2. Skewed Left Sided Difference Scheme

In skewed left sided difference scheme the number of nodes required to get fourth order approximation for a second derivative of a function is =6 ( A+O) where A is order of approximation desired and O is the order of the derivative. In this case A=4 and O=2.

Let us assume our hypothesis equation to be:

`(/_\\ x)^2 f\'\'(i)=af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)` ....(5)

We will have to expand terms of the RHS of the equation (5) that is f(i-5), f(i-4), f(i-3),f(i-2) and f(i-1) about f(i) using taylor series.

`f(i-5)=f(i)-5/_\\xf\'(i)+((5/_\\x)^2f\'\'(i))/(2!)-((5/_\\x)^3f\'\'\'(i))/(3!)+((5/_\\x)^4f\'\'\'\'(i))/(4!)-((5/_\\x)^5f\'\'\'\'\'(i))/(5!)+...`

`f(i-4)=f(i)-4/_\\xf\'(i)+((4/_\\x)^2f\'\'(i))/(2!)-((4/_\\x)^3f\'\'\'(i))/(3!)+((4/_\\x)^4f\'\'\'\'(i))/(4!)-((4/_\\x)^5f\'\'\'\'\'(i))/(5!)+...`

`f(i-3)=f(i)-3/_\\xf\'(i)+((3/_\\x)^2f\'\'(i))/(2!)-((3/_\\x)^3f\'\'\'(i))/(3!)+((3/_\\x)^4f\'\'\'\'(i))/(4!)-((3/_\\x)^5f\'\'\'\'\'(i))/(5!)...`

`f(i-2)=f(i)-2/_\\xf\'(i)+((2/_\\x)^2f\'\'(i))/(2!)-((2/_\\x)^3f\'\'\'(i))/(3!)+((2/_\\x)^4f\'\'\'\'(i))/(4!)-((2/_\\x)^5f\'\'\'\'\'(i))/(5!)...`

`f(i-1)=f(i)-/_\\xf\'(i)+((/_\\x)^2f\'\'(i))/(2!)-((/_\\x)^3f\'\'\'(i))/(3!)+((/_\\x)^4f\'\'\'\'(i))/(4!)-((/_\\x)^5f\'\'\'\'\'(i))/(5!)...`

Creating a Taylor\'s Table to evaluate the coefficients a,b,c,d,e and g for the equation (3) using the above Taylor\'s series expansions of f(i-1),f(i-2),f(i-3),f(i-4),f(i-5).

Each row of the Taylor\'s table is obtained from the coefficients of the Taylor\'s series on multiplying a with f(i-5), b with f(i-4) expansion, c with f(i-3) expansion, d with f(i-2) expansion, e with f(i-1) expansion and g with f(i).

In the Taylors table we have equated the sum of coefficients of each of the f(i), f\'(i), f\'\'\'(i), f\'\'\'\'(i), f\'\'\'\'\'(i) to zeros and f\'\'(i) to 1 because we desire to obtain an expression in terms of f\'\'(i) in the RHS of the hypothesis equation (5).

The last column of the table is not required because we need 6 equations to solve for 6 unknowns but it is shown just to show the process of filling the table.

The linear system of equations obtained from Taylor\'s table are:

a+b+c+d+e+g=0

-5a-4b-3c-2d-e=0

25a/2+8b+9c/2+2d+e/2=1

-125a/2-64b/6-27c/6-8d/6-e/6=0

625a/24+256b/24+81c/24+16d/24+e/24=0

-3125a/120-1024b/120-243c/120-32d/120-e/120=0

We can solve the system of equations using Matlab and get the value of the coefficients:

%% Computing the value of Co-ef for Skewed Left Sided Scheme
B=[0;0;1;0;0;0];
A_L=[1,1,1,1,1,1;-5,-4,-3,-2,-1,0;25/2,16/2,9/2,4/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];
X_L=A_L\\B;

The above matlab code computes for the coefficients values(a,b,c,d,e,g) in a matrix X_L. The A_L is the coefficients Matrix of the above system of linear equations, B matrix gives the residuals of the linear equations. The equation is of the form A_L*X_L=B.

So, X_L= (Inversion of A_L matrix) * B

The output from Matlab is:

So the hypothesis equation for LDS scheme can be approximately written as:

`(/_\\ x)^2 f\'\'(i)=-(5f(i-5))/6+(61f(i-4))/12-13f(i-3)+(107f(i-2))/6-(77f(i-1))/6+(15f(i))/4+O(/_\\x)^6`

`f\'\'(i)=(-(5f(i-5))/6+(61f(i-4))/12-13f(i-3)+(107f(i-2))/6-(77f(i-1))/6+(15f(i))/4)/(/_\\x)^2+O(/_\\x)^4`

`f\'\'(i)~~(45f(i)-154f(i-1)+214f(i-2)-156f(i-3)+61f(i-4)-10f(i-5))/(12(/_\\x)^2)`... (6)

We will use this equation (6) in Matlab program to compute Skewed Left sided Difference Scheme fourth order approximation of second derivative and corresponding error from exact second derivative.

Function code for Central Difference Scheme

function result= central_diff_fourth_order(x,dx)
% Function is f(x)=e^x*cosx
result=(-exp(x-2*dx)*cos(x-2*dx)+16*exp(x-dx)*cos(x-dx)-30*exp(x)*cos(x)+16*exp(x+dx)*cos(x+dx)+...
        -exp(x+2*dx)*cos(x+2*dx))/(12*dx^2);
end
% f\'\'(x)=(-f(i-2)+16f(i-1)-30f(i)+16f(i+1)-f(i+2))/12dx^2

Function code for Skewed Right sided Difference Scheme

function result= right_skewed_fourth_order(x,dx)
% Function is f(x)=e^x*cosx
result=(45*exp(x)*cos(x)-154*exp(x+dx)*cos(x+dx)+214*exp(x+2*dx)*cos(x+2*dx)-156*exp(x+3*dx)*cos(x+3*dx)+...
        +61*exp(x+4*dx)*cos(x+4*dx)-10*exp(x+5*dx)*cos(x+5*dx))/(12*dx^2);
end
% f\'\'(x)=(45f(i)-154f(i+1)+214f(i+2)-156f(i+3)+61f(i+4)-10f(i+5))/12dx^2
    

Function code for Skewed Left sided Difference Scheme

function result= left_skewed_fourth_order(x,dx)
% Function is f(x)=e^x*cosx
result=(45*exp(x)*cos(x)-154*exp(x-dx)*cos(x-dx)+214*exp(x-2*dx)*cos(x-2*dx)-156*exp(x-3*dx)*cos(x-3*dx)+...
        +61*exp(x-4*dx)*cos(x-4*dx)-10*exp(x-5*dx)*cos(x-5*dx))/(12*dx^2);
end
% f\'\'(x)=(-10f(i-5)+61f(i-4)-156f(i-3)+214f(i-2)-154f(i-1)+45f(i))/12dx^2
    

Main Program

%% Computing fourth Order Approximation of Second Derivative for f(x)=e^x*cos(x)
% Variables
clc
close all
clear all
x=pi/3;
dx=linspace(pi/20000,pi/10,200);
Exact_diff= -2*exp(x)*sin(x);

%% Computing Numerically using discretization

for i=1 :length(dx)

 Central_diff(i)=central_diff_fourth_order(x,dx(i));
 
 Right_skewed(i)=right_skewed_fourth_order(x,dx(i));
 
 Left_skewed(i)=left_skewed_fourth_order(x,dx(i));
 
end
 
 %% Error Calculation
 
 Err_central=abs(Exact_diff-Central_diff);
 Err_right=abs(Exact_diff-Right_skewed);
 Err_left= abs(Exact_diff-Left_skewed);
 
 %% PLOTS
 
 figure(1);
 plot(dx,Err_central,dx,Err_right,dx,Err_left,\'linewidth\',2);
 legend({\'Central Diff\',\'Skewed Right sided Diff\',\'Skewed Left sided Diff\'},\'Location\',\'northwest\');
 xlabel(\'Step-size (dx)\');
 ylabel(\'Error\');
 title(\'Plot of dx Vs Error of Different 4th order Approximation Schemes\');
 figure(2);
 loglog(dx,Err_central,dx,Err_right,dx,Err_left,\'linewidth\',2);
 legend({\'Central Diff\',\'Skewed Right sided Diff\',\'Skewed Left sided Diff\'},\'Location\',\'northwest\');
 xlabel(\'Step-size (dx)\');
 ylabel(\'Error\');
 title(\'Plot of dx Vs Error of Different 4th order Approximation Schemes\');

PLOTS

From the plot above we can see that the Central Difference Scheme shows least error and seems flat followed by Skewed Left sided sheme and then Right sided scheme which shows a sharp curved decrease. But for lower values of dx we cannot make any inference so we will plot again in logarithm scale using log plots.

From the plots above we can say that again Central Difference Scheme shows least error compared to Skewed schemes. Between dx values of the approximately of the order of 10^-2 to 10^-1 the Left and Right sided Skemes behave similarly and overlaps wih each other.

Another important conclusion that can be made is that in all the three schemes error start to increase as dx becomes very small closer to 10^-3. This anamoly behaviour is due to machine precision(matlab precision) which tends to round off the error when its very small.

For values of dx lower than 10^-3 the left sided scheme shows less error than right sided sheme. This is also true for values of dx closer or greater than 10^-1.

Skewed schemes usefulness

Skewed shemes are useful most of the time as it doesnot require information from its both sides of the node ( that too symmetric position) about which we tend to compute any variable which is the primary condition for central shemes. So if we have information on one side(either left or right only) of the node then we can compute the desired parameter at that node using information from one sided nodes.

Consider the best examples where we can use skewed schemes are hyperbolic (say invicid supersonic flow) or parabolic(viscous flow) partial differential equations. The solutions of these equations are evaluated by a marching procedure which require information from previous nodes and always move in forward direction in time or space to compute.