Numerical comparison of the absolute error of various 4th order approximations of the 2nd order derivative of a function
Part1:
Here we shall try to derive the 4th order approximations of the 2nd order derivative by using following methods:
- Central difference
- Skewed right sided difference
- Skewed left sided difference
For the above analysis we will use the following table which shows the coefficients of taylor-series expansion.
Table 1:
| f(i-5) | f(i-4) | f(i-3) | f(i-2) | f(i-1) | f(i) | f(i+1) | f(i+2) | f(i+3) | f(i+4) | f(i+5) | |
| f(i) | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| f1(i)*∆x | -5 | -4 | -3 | -2 | -1 | 0 | 1 | 2 | 3 | 4 | 5 |
| f2(i)*∆x^2 | 5^2/2 | 4^2/2 | 3^2/2 | 2^2/2 | 1/2 | 0 | 1/2 | 2^2/2 | 3^2/2 | 4^2/2 | 5^2/2 |
| f3(i)*∆x^3 | -5^3/6 | -4^3/6 | -3^3/6 | -2^3/6 | -1/6 | 0 | 1/6 | 2^3/6 | 3^3/6 | 4^3/6 | 5^3/6 |
| f4(i)*∆x^4 | 5^4/24 | 4^4/24 | 3^4/24 | 2^4/24 | 1/24 | 0 | 1/24 | 2^4/24 | 3^4/24 | 4^4/24 | 5^4/24 |
| f5(i)*∆x^5 | -5^5/120 | -4^5/120 | -3^5/120 | -2^5/120 | -1/120 | 0 | 1/120 | 2^5/120 | 3^5/120 | 4^5/120 | 5^5/120 |
| f6(i)*∆x^6 | -5^6/720 | -4^6/720 | -3^6/720 | -2^6/720 | -1/720 | 1 | 1/720 | 2^5/720 | 3^5/720 | 4^5/720 | 5^5/720 |
1.Central differencing Method:
Lets us try to approcimate, `(del^2f)/(delx^2)` using polynomials, as follows:
`(del^2f)/(delx^2) approx af(i-2) + bf(i-1) + cf(i) + df(i+1) + ef(i+2)`
where,
We can reduce Table 1 to the following:
| `f(i-2)` | `f(i-1)` | `f(i)` | `f(i+1)` | `f(i+2)` | |
| `f(i)` | 1 | 1 | 1 | 1 | 1 |
| `f\'(i)∆x` | -2 | -1 | 0 | 1 | 2 |
| `f\'\'(i)∆x^2` | 2^2/2 | 1/2 | 0 | 1/2 | 2^2/2 |
| `f\'\'\'(i)∆x^3` | -2^3/6 | -1/6 | 0 | 1/6 | 2^3/6 |
| `f^(iv)(i)∆x^4` | 2^4/24 | 1/24 | 0 | 1/24 | 2^4/24 |
| `f^(v)(i)∆x^5` | -2^5/120 | -1/120 | 0 | 1/120 | 2^5/120 |
so to obtain 2nd order derivative, we will form a equation, AX=B, where;
A=
| 1 | 1 | 1 | 1 | 1 |
| -2 | -1 | 0 | 1 | 2 |
| 2^2/2 | 1/2 | 0 | 1/2 | 2^2/2 |
| -2^3/6 | -1/6 | 0 | 1/6 | 2^3/6 |
| 2^4/24 | 1/24 | 0 | 1/24 | 2^4/24 |
X=
| a |
| b |
| c |
| d |
| e |
B=
| 0 |
| 0 |
| 1 |
| 0 |
| 0 |
upon solving and rationalising we get,
| a | -0 + 1/(-12) |
| b | 1 + 1/(3) |
| c | -3 + 1/(2) |
| d | 1 + 1/(3) |
| e | -0 + 1/(-12) |
Hence, negleting h.o.t we can write the 4th order central difference approxiamtion of 2nd order derivative as follows;
`\\frac{\\partial^{(2)}f}{\\partial x^{(2)}}\\approx\\frac{-f(x-2Deltax)+16f(x-1Deltax)-30f(x)+16f(x+Deltax)-f(x+2Deltax)}{12Deltax^{2}}`
2.Skewed right sided difference
Lets us try to approcimate, `(del^2f)/(delx^2)` using polynomials, as follows:
`(del^2f)/(delx^2) approx af(i) + bf(i+1) + cf(i+2)+df(i+3) + ef(i+4) + gf(i+5)`
where,
We can reduce Table 1 to the following:
| f(i) | f(i+1) | f(i+2) | f(i+3) | f(i+4) | f(i+5) | |
| f(i) | 1 | 1 | 1 | 1 | 1 | 1 |
| f1(i)*∆x | 0 | 1 | 2 | 3 | 4 | 5 |
| f2(i)*∆x^2 | 0 | 1/2 | 2^2/2 | 3^2/2 | 4^2/2 | 5^2/2 |
| f3(i)*∆x^3 | 0 | 1/6 | 2^3/6 | 3^3/6 | 4^3/6 | 5^3/6 |
| f4(i)*∆x^4 | 0 | 1/24 | 2^4/24 | 3^4/24 | 4^4/24 | 5^4/24 |
| f5(i)*∆x^5 | 0 | 1/120 | 2^5/120 | 3^5/120 | 4^5/120 | 5^5/120 |
so to obtain 2nd order derivative, we will form a equation, AX=B, where;
A=
| 1 | 1 | 1 | 1 | 1 | 1 |
| 0 | 1 | 2 | 3 | 4 | 5 |
| 0 | 1/2 | 2^2/2 | 3^2/2 | 4^2/2 | 5^2/2 |
| 0 | 1/6 | 2^3/6 | 3^3/6 | 4^3/6 | 5^3/6 |
| 0 | 1/24 | 2^4/24 | 3^4/24 | 4^4/24 | 5^4/24 |
| 0 | 1/120 | 2^5/120 | 3^5/120 | 4^5/120 | 5^5/120 |
X=
| a |
| b |
| c |
| d |
| e |
| g |
B=
| 0 |
| 0 |
| 1 |
| 0 |
| 0 |
| 0 |
upon solving and rationalising we get,
| a | -1 + 1/(6) |
| b | 5 + 1/(12) |
| c | -13 |
| d | 18 + 1/(-6) |
| e | -13 + 1/(6) |
| g | 4 + 1/(-4) |
Hence, negleting h.o.t we can write the 4th order central difference approxiamtion of 2nd order derivative as follows;
`\\frac{\\partial^{(2)}f}{\\partial x^{(2)}}\\approx\\frac{45f(x)-154f(x+Deltax)+214f(x+2Deltax)-156f(x+3Deltax)+61f(x+4Deltax)-10f(x+5Deltax)}{12Deltax^{2}}`
3.Skewed left sided difference
Lets us try to approcimate, `(del^2f)/(delx^2)` using polynomials, as follows:
`(del^2f)/(delx^2) approx af(i-5) + bf(i-4) + cf(i-3)+df(i-2) + ef(i-1) + gf(i)`
where,
We can reduce Table 1 to the following:
| f(i) | f(i+1) | f(i+2) | f(i+3) | f(i+4) | f(i+5) | |
| f(i) | 1 | 1 | 1 | 1 | 1 | 1 |
| f1(i)*∆x | -5 | -4 | -3 | -2 | -1 | 0 |
| f2(i)*∆x^2 | 5^2/2 | 4^2/2 | 3^2/2 | 2^2/2 | 1/2 | 0 |
| f3(i)*∆x^3 | -5^3/6 | -4^3/6 | -3^3/6 | -2^3/6 | -1/6 | 0 |
| f4(i)*∆x^4 | 5^4/24 | 4^4/24 | 3^4/24 | 2^4/24 | 1/24 | 0 |
| f5(i)*∆x^5 | -5^5/120 | -4^5/120 | -3^5/120 | -2^5/120 | -1/120 | 0 |
so to obtain 2nd order derivative, we will form a equation, AX=B, where;
A=
| 1 | 1 | 1 | 1 | 1 | 1 |
| -5 | -4 | -3 | -2 | -1 | 0 |
| 5^2/2 | 4^2/2 | 3^2/2 | 2^2/2 | 1/2 | 0 |
| -5^3/6 | -4^3/6 | -3^3/6 | -2^3/6 | -1/6 | 0 |
| 5^4/24 | 4^4/24 | 3^4/24 | 2^4/24 | 1/24 | 0 |
| -5^5/120 | -4^5/120 | -3^5/120 | -2^5/120 | -1/120 | 0 |
X=
| a |
| b |
| c |
| d |
| e |
| g |
and B=
| 0 |
| 0 |
| 1 |
| 0 |
| 0 |
| 0 |
upon solving and rationalising we get,
| a | -1 + 1/(6) |
| b | 5 + 1/(12) |
| c | -13 |
| d | 18 + 1/(-6) |
| e | -13 + 1/(6) |
| g | 4 + 1/(-4) |
Hence, negleting h.o.t we can write the 4th order central difference approxiamtion of 2nd order derivative as follows;
`\\frac{\\partial^{(2)}f}{\\partial x^{(2)}}\\approx\\frac{-10f(x-5Deltax)+61f(x-4Deltax)-156f(x-3Deltax)+214f(x-2Deltax)-154f(x-Deltax)+45f(x)}{12Deltax^{2}}`
Part2:
Matlab program to evaluate the second order derivative of the function `e^(x)cos(x)`analytically and approximately with the 3 numerical approximations(4th order), which was derived previously.
Functions:
1.central_scheme.m
function cd_out=central_scheme(x,dx)
%central_differencing--------------------
%d^2y/dx^2= (-1*f[i-2]+16*f[i-1]-30*f[i]+16*f[i+1]-1*f[i+2])/(12*dx^2)
%----------------------------------------
cd_out=(-(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));
end2. skewed_right_scheme.m
function sr_out=skewed_right_scheme(x,dx)
%skewed_right_differencing---------------
%d^2y/dx^2= (45*f[i]-154*f[i+1]+214*f[i+2]-156*f[i+3]+61*f[i+4]-10*f[i+5])/(12*dx^2)
sr_out= ((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));
end3. skewed_left_scheme.m
function sl_out=skewed_left_scheme(x,dx)
%skewed_left_differencing----------------
%d^2y/dx^2 = (-10*f[i-5]+61*f[i-4]-156*f[i-3]+214*f[i-2]-154*f[i-1]+45*f[i+0])/(12*dx^2)
sl_out=((-10*(exp(x-(5*dx))*cos(x-(5*dx))))+ (61*(exp(x-(4*dx))*cos(x-(4*dx))))- (156*(exp(x-(3*dx))*cos(x-(3*dx))))+ (214*(exp(x-(2*dx))*cos(x-(2*dx))))-(154*(exp(x-(dx))*cos(x-(dx))))+(45*(exp(x)*cos(x))))/(12*(dx^2));
end4. fd_scheme.m
function fd_out=fd_scheme(x,dx)
%forward_differencing----------------
%d^2y/dx^2=(1*f[i+0]-2*f[i+1]+1*f[i+2])/(dx^2)
fd_out=((exp(x)*cos(x))-(2*(exp(x+(dx))*cos(x+(dx)))) + (exp(x+(2*dx))*cos(x+(2*dx))))/(dx^2);
end5. bd_scheme.m
function bd_out=bd_scheme(x,dx)
%forward_differencing----------------
%d^2y/dx^2=(1*f[i-2]-2*f[i-1]+1*f[i+0])/(dx^2)
bd_out=((exp(x-(2*dx))*cos(x-(2*dx)))-(2*(exp(x-(dx))*cos(x-(dx))))+((exp(x)*cos(x))))/(dx^2);
endScript:
%-----Differencing Methods-----
%Author: Mahammad Sehzad Alli
%---------------------------------
%Given function is (e^x)cos(x)
%-----------------------------------------
clc
close all % deletes all figures whose handles are not hidden.
clear % clear removes all variables from the current workspace
%--------Input----------------------------
x=linspace(pi/2,2*pi,20);
dx=[pi/6,pi/10,pi/20];
ad=zeros(numel(x),3);
cd=zeros(numel(x),3);
sr=zeros(numel(x),3);
sl=zeros(numel(x),3);
%-----------------------------------------
%analytical_derivative-------------------
for i=1:numel(x)
for n=1:numel(dx)
ad(i,n)= -2*exp(x(i))*sin(x(i));
cd(i,n)=central_scheme(x(i),dx(n));
sr(i,n)=skewed_right_scheme(x(i),dx(n));
sl(i,n)=skewed_left_scheme(x(i),dx(n));
end
end
%---------------------
figure(1)
plot(x,ad(:,2),\'-s\',x,cd(:,1),\'-*\',x,sr(:,1),\'-+\',x,sl(:,1),\'-x\',\'LineWidth\',1)
xlabel \'x\'
ylabel (\'d^{2}y/dx^{2}\')
legend(\'Analytical\',\'Central\',\'Skewed right\',\'Skewed left\',\'location\',\'southwest\')
title(\'Fourth Order Approximation\')
axis([1.5,6.5,-800,400])
txt = [\'dx=\',num2str(dx(1))];
text(4,-600,txt)
%-------------Saving------------------------------------
saveas(1,\'solution_comparision_1.png\')
%---------------------
figure(2)
plot(x,ad(:,2),\'-s\',x,cd(:,2),\'-*\',x,sr(:,2),\'-+\',x,sl(:,2),\'-x\',\'LineWidth\',1)
xlabel \'x\'
ylabel (\'d^{2}y/dx^{2}\')
legend(\'Analytical\',\'Central\',\'Skewed right\',\'Skewed left\',\'location\',\'southwest\')
title(\'Fourth Order Approximation\')
axis([1.5,6.5,-800,400])
txt = [\'dx=\',num2str(dx(2))]; %else str=sprintf(\'dx=%f\',dx(1))
text(4,-600,txt)
%-------------Saving------------------------------------
saveas(2,\'solution_comparision_2.png\')
%---------------------
figure(3)
plot(x,ad(:,3),\'-s\',x,cd(:,3),\'-*\',x,sr(:,3),\'-+\',x,sl(:,3),\'-x\',\'LineWidth\',1)
xlabel \'x\'
ylabel (\'d^{2}y/dx^{2}\')
legend(\'Analytical\',\'Central\',\'Skewed right\',\'Skewed left\',\'location\',\'southwest\')
title(\'Fourth Order Approximation\')
axis([1.5,6.5,-800,400])
txt = [\'dx=\',num2str(dx(3))];
text(4,-600,txt)
%-------------Saving------------------------------------
saveas(3,\'solution_comparision_3.png\')Output:
Comment:
1.As the spacing between nodes (dx) decreases solution accuracy increases.
2. At higher dx, accuracy of skewed right sided difference scheme decreases drastically.
Matlab program to evaluate and plot the comparision of the absolute error between the above mentioned schemes(3 numerical approximations of 4th order) and analytical derivation of second order derivative of the function (`e^(x)cos(x)`)
Function
The above mentioned functions used in this script too.
6.differencing_error.m
function [c,sr,sl,fd,bd]=differencing_error(x,dx)
%-----------------------------------------
analytical_derivative= -2*exp(x)*sin(x);
%--------/Fourth order approximation\\--------
%---1.Central Diference-----------------------------
central_differencing=central_scheme(x,dx);
%-------error-----------
c=abs(central_differencing-analytical_derivative);
%---2.Skewed right sided difference-----------------
skewed_right_differencing=skewed_right_scheme(x,dx);
%-------error-----------
sr=abs(skewed_right_differencing-analytical_derivative);
%---3.Skewed left sided difference-----------------
skewed_left_differencing=skewed_left_scheme(x,dx);
%-------error------------
sl=abs(skewed_left_differencing-analytical_derivative);
%---4.forward difference-----------------
forward_differencing=fd_scheme(x,dx);
%-------error------------
fd=abs(forward_differencing-analytical_derivative);
%---5.backward difference-----------------
backward_differencing=bd_scheme(x,dx);
%-------error------------
bd=abs(backward_differencing-analytical_derivative);
Script
%-----Error Comparision-----
%Author: Mahammad Sehzad Alli
%---------------------------------
%Given function is (e^x)cos(x)
%-----------------------------------------
clc
close all % deletes all figures whose handles are not hidden.
clear % clear removes all variables from the current workspace
%--------Input----------------------------
x=pi/2;
dx=linspace(pi/3,pi/40,200)\';
er_dif=zeros(numel(dx):1,3);
%-----------------------------------------
for i=1:numel(dx)
[er_dif(i,1),er_dif(i,2),er_dif(i,3),er_dif(i,4),er_dif(i,5)]=differencing_error(x,dx(i));
end
%-------------Figure---------------------------------
figure(1)
loglog(dx,er_dif(:,1),dx,er_dif(:,2),dx,er_dif(:,3),dx,er_dif(:,4),\'--\',dx,er_dif(:,5),\'-\',\'LineWidth\',1.5)
xlabel \'dx\'
ylabel \'Absolute error magnitude\'
legend(\'Central\',\'Skewed right\',\'Skewed left\',\'FDS\',\'BDS\',\'location\',\'Northwest\')
title(\'Error comparision\')
saveas(1,\'error_comparision_1.png\')
Output:
Comment:
1. Central difference scheme offers better accuracy compared to skewed schemes.
2. F.D.S & B.D.S scheme are first order schemes, hence their accuracy is less compared to all fourth order schemes.
3. At very less dx, value round-off error becomes significant, quite apparant from increase in the absolute error value of both the skewed difference schemes.
Part3
why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?
Let say we are trying to solve a Hagen–Poiseuille flow,(Laminar flow of an incompressible and Newtonian fluid, through a long cylindrical pipe of constant cross section).
If we try to find the flow field information (let say Temperature, T) at point B, using central approximation, we may have sufficient information available at the neighbouring region to compute it.
But, at Point B and C, we have disadvantage of lack of sufficient information in adjacent region.
However, we can overcome this by using skewed scheme,
In certain fluid flow problems, like advection, due to lack of flow data at previous position( point A) or later position (point C) it is imperative to use the flow data available upstream or downstream to compute the flow data of current node. Hence this methodology is called as skewed finite difference scheme.