Derive the 4th order approximations of the second order derivative using various finite difference schemes .
In this project we are going to derive the 4th order approximations of the second-order derivative using the following schemes :
1.) Central Differencing
2.) Skewed right-sided difference
3.) Skewed left-sided difference
The analytical function for which we will derive these schemes is f(x) = `exp(x)*cos(x)`.
An insight on Central differencing: The central differencing scheme is a finite difference method employed in calculating the approximations of the derivatives of the analytical functions by using information from both sides of the point. Consider the figure shown here :
In the central differencing scheme, the information is taken from both sides of the point.
An insight on Skewed right-sided differencing: The skewed right-sided differencing scheme ( forward differencing ) is a finite difference method employed in calculating the approximations of the derivatives of the analytical functions by using information from the right side of the point only. Consider the figure shown here :
In the skewed right-sided differencing scheme the information is taken from the right side of the point.
An insight on Skewed left-sided differencing: The skewed left-sided differencing scheme ( backward differencing ) is a finite difference method employed in calculating the approximations of the derivatives of the analytical functions by using information from the left side of the point only. Consider the figure shown here :
In the skewed left-sided differencing scheme the information is taken from the left side of the point.
To derive the information for the 4th order approximation of the second-order derivative, the Taylor table method is employed.
Consider the derivatives are to be evaluated at a point x and the discrete steps available are dx.
Central Differencing
As in the central differencing, the information is to be taken from both the sides of the point equally, consider the following Taylor table :
Now the second derivative with central differencing scheme can be represented as :
`f''(x) = af(x-2dx)+bf(x-dx)+cf(x)+df(x+dx)+ef(x+2dx)`
So, multiplying the appropriate co-efficients in the TAYLOR TABLE we get :
On totalling all the entries of the table column wise we get :
1.) `f''(x) xx (dx)^2= af(x-2dx)+bf(x-dx)+cf(x)+df(x+dx)+ef(x+2dx)`
2.) `a+b+c+d+e = 0`
3.) `-2a - b+0+d+2e = 0`
4.) `2a + frac(b)(2) + 0 +frac(d)(2) + 2e = 1`
5.) `frac(-8a)(3) + frac(-b)(6) + 0+frac(d)(6)+frac(8e)(3) = 0`
6.) `frac(2a)(3) + frac(b)(24) +0 +d/24 + frac(2e)(3) = 0`
Solving these equations we will get the values of a, b, c, d, and e. MATLAB SYMS TOOLBOX is used to solve these equations as follows :
clear all
close all
clc
% I will represent a matrix of unknown variables using syms toolbox
syms a b c d e;
% The equations are as follows:
% a + b + c + d + e = 0 ;
% -2*a - b + 0 + d + 2*e = 0;
% 2a + b/2 + 0 + d/2 + 2e = 1/dx^2;
% -8a/3 + -b/6 + 0 + d/6 + 8e/3 = 0;
% 2a/3 + b/24 + 0 + d/24 + 2e/3 = 0;
% Representing the equations in the form of a matrix AX = B
A = [ 1 1 1 1 1 ; ...
-2 -1 0 1 2 ; ...
2 1/2 0 1/2 2 ; ...
-8/3 -1/6 0 1/6 8/3 ; ...
2/3 1/24 0 1/24 2/3 ];
X = [ a ; ...
b ; ...
c ; ...
d ; ...
e ] ;
B = [ 0 ; ...
0 ; ...
1 ; ...
0 ; ...
0];
% Solving for x we do X = inverse(A)*B
C = inv(A)*B ;
a = C(1) ;
b = C(2) ;
c = C(3) ;
d = C(4) ;
e = C(5) ;
The coefficients a, b, c, d, e are determined for the forward differencing scheme as shown above.
a = -0.0833
b = 1.3333
c = -2.5000
d = 1.3333
e = -0.0833
Therefore we have :
`f''(x) = frac(af(x-2dx)+bf(x-dx)+cf(x)+df(x+dx)+ef(x+2dx))(dx^2)`
Switching to MATLAB and writing the code to compute this derivative :
% The given analytical function is f(x) = exp(x)*cos(x)
% The analytical derivative is given by
%
% f'(x) = exp(x)cos(x)-exp(x)sin(x)
%
% I will compute the derivative at x = pi/3 ;
% Defining the x-value at which the derivative is to be calculated
x = pi/3;
% Defining the dx parameter as a matrix so that we can see compare the
% error value with that of dx
dx = [ pi/3 ; pi/30 ; pi/300 ; pi/3000 ; pi/30000];
% The analytical derivative first order
analytical_derivative_first_order = exp(x)*cos(x) - exp(x)*sin(x);
% The analytical derivative second order
% analytical_derivative_second_order = exp(x)*cos(x) - exp(x)*sin(x) - exp(x)*sin(x) - exp(x)*cos(x) ;
analytical_derivative_second_order = -2*exp(x)*sin(x);
% Computing the numerical second order derivative of fourth order
% Approximation using the central differencing scheme just derived
numerical_central_scheme_1 = a*exp(x-2*dx).*cos(x-2*dx) + b*exp(x-dx).*cos(x-dx)...
+ c*exp(x).*cos(x) + d*exp(x+dx).*cos(x+dx)...
+ e*exp(x+2*dx).*cos(x+2*dx);
numerical_derivative_central_scheme = numerical_central_scheme_1./(dx).^2;
% Error in the numerical vs. analytical derivative
error_central_scheme = abs(numerical_derivative_central_scheme-analytical_derivative_second_order);
% plotting of the error
plot(dx,error_central_scheme,'ks-','markersize',12,'markerfacecolor','y');
xlim([-0.1 1.2]);
ylim([-0.01 0.33]);
xlabel(' dx values ');
ylabel(' Error');
title(' Error in central differencing scheme');
The graph representing the errors is shown below :
Forward Differencing Scheme
As in the forward differencing ( skewed right-sided ), the information is to be taken from the right side of the point only, consider the following Taylor table :
Now the second derivative with forward differencing scheme ( right skewed differencing scheme ) can be represented as :
`f''(x) = af(x)+bf(x+dx)+cf(x+2dx)+df(x+3dx)+ef(x+4dx)`
So, multiplying the appropriate coefficients in the TAYLOR TABLE we get :
On totalling all the entries of the table column wise we get :
1.) `af(x) +bf(x+dx)+cf(x+2dx)+df(x+3dx) + ef(x+4dx) = f''(x)xx (dx)^2`
2.) `b + 2c + 3d +4e =0`
3.) `a+b+c+d+e=0`
4.) `b/2 + 2c + frac(9d)(2) +8e = 1`
5.) `b/6 + frac(4c)(3) + frac(9d)(2)+frac(32e)(3) = 0`
6.) `b/24 + frac(2c)(3) +frac(27d)(8)+frac(32e)(3) =0`
Solving these equations we will get the values of a, b, c, d, and e. MATLAB SYMS TOOLBOX is used to solve these equations as follows :
clear all
close all
clc
% I will represent a matrix of unknown variables using syms toolbox
syms a b c d e;
% the equations are as follows:
% a + b + c + d + e = 0 ;
% 0 + b + 2c + 3d + 4e = 0;
% 0 + b/2 + 2c + 9d/2 +8e = 1;
% 0 + b/6 + 8c/3 + 9d/2 + 32e/3 = 0;
% 0 + b/24 + 2c/3 + 27d/8 + 32e/3 = 0;
% representing the equations in the form of a matrix AX = B
A = [ 1 1 1 1 1 ; ...
0 1 2 3 4 ; ...
0 1/2 2 9/2 8 ; ...
0 1/6 4/3 9/2 32/3 ; ...
0 1/24 2/3 27/8 32/3 ];
X = [ a ; ...
b ; ...
c ; ...
d ; ...
e ] ;
B = [ 0 ; ...
0 ; ...
1 ; ...
0 ; ...
0];
% for solving for x we do X = inverse(A)*B
C = inv(A)*B ;
a = C(1) ;
b = C(2) ;
c = C(3) ;
d = C(4) ;
e = C(5) ;
The coefficients a, b, c, d, e are determined for the forward differencing scheme as shown above.
a = 2.9167
b = -8.667
c = 9.5000
d = -4.6667
e = 0.9167
Therefore we have :
`f''(x) = frac(af(x)+bf(x+dx)+cf(x+2dx)+df(x+3dx)+ef(x+4dx))(dx^2)`
Switching to MATLAB and writing the code to compute this derivative :
% the given analytical function is f(x) = exp(x)*cos(x)
% The analytical derivative is given by
%
% f'(x) = exp(x)cos(x)-exp(x)sin(x)
%
% I will compute the derivative at x = pi/3 ;
% defining the x-value at which the derivative is to be calculated
x = pi/3;
% Defining the dx parameter as a matrix so that we can see compare the
% error value with that of dx
dx = [ pi/3 ; pi/30 ; pi/300 ; pi/3000 ; pi/30000];
% The analytical derivative first order
analytical_derivative_first_order = exp(x)*cos(x) - exp(x)*sin(x);
% The analytical derivative second order
% analytical_derivative_second_order = exp(x)*cos(x) - exp(x)*sin(x) - exp(x)*sin(x) - exp(x)*cos(x) ;
analytical_derivative_second_order = -2*exp(x)*sin(x);
% computing the numerical second order derivative of fourth order
% approximation using the central differencing scheme just derived
numerical_forward_differencing_scheme_1 = a*exp(x).*cos(x) + b*exp(x+dx).*cos(x+dx)...
+ c*exp(x+2*dx).*cos(x+2*dx) + d*exp(x+3*dx).*cos(x+3*dx)...
+ e*exp(x+4*dx).*cos(x+4*dx);
numerical_derivative_forward_scheme = numerical_forward_differencing_scheme_1./(dx).^2;
% Error in the numerical vs. analytical derivative
error_forward_scheme = abs(numerical_derivative_forward_scheme-analytical_derivative_second_order);
% plotting of the error
plot(dx,error_forward_scheme,'ks-','markersize',12,'markerfacecolor','y');
xlim([-0.1 1.2]);
ylim([-1 61]);
xlabel(' dx values ');
ylabel(' Error');
title(' Error in central differencing scheme');
The graph representing the errors in the forward differencing is shown below :
Backward Differencing Scheme
As in the backward differencing ( skewed left-sided ) the information is to be taken from the left side of the point only, consider the following taylor table :
Now the second derivative with forward differencing scheme ( right skewed differencing scheme ) can be represented as :
`f''(x) = af(x)+bf(x-dx)+cf(x-2dx)+df(x-3dx)+ef(x-4dx)`
So, multiplying the appropriate coefficients in the TAYLOR TABLE we get :
On totalling all the entries of the table column wise we get :
1.) `af(x) +bf(x-dx)+cf(x-2dx)+df(x-3dx) + ef(x-4dx) = f''(x)xx(dx)^2`
2.) `b + 2c + 3d +4e =0`
3.) `a+b+c+d+e=0`
4.) `b/2 + 2c + frac(9d)(2) +8e = 1`
5.) `b/6 + frac(4c)(3) + frac(9d)(2)+frac(32e)(3) = 0`
6.) `b/24 + frac(2c)(3) +frac(27d)(8)+frac(32e)(3) =0`
Solving these equations we will get the values of a, b, c, d, and e. MATLAB SYMS TOOLBOX is used to solve these equations as follows :
clear all
close all
clc
% I will represent a matrix of unknown variables using syms toolbox
syms a b c d e;
% the equations are as follows:
% a + b + c + d + e = 0 ;
% 0 + b + 2c + 3d + 4e = 0;
% 0 + b/2 + 2c + 9d/2 +8e = 1;
% 0 + b/6 + 8c/3 + 9d/2 + 32e/3 = 0;
% 0 + b/24 + 2c/3 + 27d/8 + 32e/3 = 0;
% representing the equations in the form of a matrix AX = B
A = [ 1 1 1 1 1 ; ...
0 1 2 3 4 ; ...
0 1/2 2 9/2 8 ; ...
0 1/6 4/3 9/2 32/3 ; ...
0 1/24 2/3 27/8 32/3 ];
X = [ a ; ...
b ; ...
c ; ...
d ; ...
e ] ;
B = [ 0 ; ...
0 ; ...
1 ; ...
0 ; ...
0];
% for solving for x we do X = inverse(A)*B
C = inv(A)*B ;
a = C(1) ;
b = C(2) ;
c = C(3) ;
d = C(4) ;
e = C(5) ;
The coefficients a, b, c, d, e are determined for the forward differencing scheme as shown above.
a = 2.9167
b = -8.667
c = 9.5000
d = -4.6667
e = 0.9167
Therefore we have :
`f''(x) = frac(af(x)+bf(x+dx)+cf(x+2dx)+df(x+3dx)+ef(x+4dx))(dx^2)`
Switching to MATLAB and writing the code to compute this derivative :
% the given analytical function is f(x) = exp(x)*cos(x)
% The analytical derivative is given by
%
% f'(x) = exp(x)cos(x)-exp(x)sin(x)
%
% I will compute the derivative at x = pi/3 ;
% defining the x-value at which the derivative is to be calculated
x = pi/3;
% Defining the dx parameter as a matrix so that we can see compare the
% error value with that of dx
dx = [ pi/3 ; pi/30 ; pi/300 ; pi/3000 ; pi/30000];
% The analytical derivative first order
analytical_derivative_first_order = exp(x)*cos(x) - exp(x)*sin(x);
% The analytical derivative second order
% analytical_derivative_second_order = exp(x)*cos(x) - exp(x)*sin(x) - exp(x)*sin(x) - exp(x)*cos(x) ;
analytical_derivative_second_order = -2*exp(x)*sin(x);
% computing the numerical second order derivative of fourth order
% approximation using the central differencing scheme just derived
numerical_backward_differencing_scheme_1 = a*exp(x).*cos(x) + b*exp(x+dx).*cos(x+dx)...
+ c*exp(x+2*dx).*cos(x+2*dx) + d*exp(x+3*dx).*cos(x+3*dx)...
+ e*exp(x+4*dx).*cos(x+4*dx);
numerical_derivative_backward_scheme = numerical_backward_differencing_scheme_1./(dx).^2;
% Error in the numerical vs. analytical derivative
error_backward_scheme = abs(numerical_derivative_backward_scheme-analytical_derivative_second_order);
% plotting of the error
plot(dx,error_backward_scheme,'ks-','markersize',12,'markerfacecolor','y');
xlim([-0.1 1.2]);
ylim([-1 61]);
xlabel(' dx values ');
ylabel(' Error');
title(' Error in central differencing scheme');
The graph representing the errors in backward differencing is shown below :
Conclusion :
It can be observed that the error in the forward and the backward differencing schemes are quite hight for a smaller value of the discretized value of dx as compared to central differencing. Therefore if the computation power of a computer is not so good then one should proceed with a smaller value of dx and use central differencing to obtain accurate results. Sometimes it is possible that at the boundary of a mesh we need to find a value using any difference scheme, in here the forward or the backward difference scheme comes in handy as the data at the boundary point is not available in both directions i.e. it is available in one direction only.
Also, the central difference is more accurate as compared to FDS and BDS because in central difference the information is obtained from the points in the vicinity of the given point and hence the results are much more accurate.
Discussion :
The finite difference method is very useful in modern-day engineering applications like computational fluid dynamics, Structural Analysis, etc and it can solve hard engineering problems up to good accuracy which can't be solved conventionally. The method to apply the finite difference method is a bit tricky as sometimes compromising the computational power of the computer we get not so accurate results by using some other difference scheme or vice versa. It is possible that the simulation can give two entirely different results when two different types of differencing schemes are applied, so we must choose carefully in applying the differencing schemes.