Analysis of fourth order accurate finite difference schemes using Taylor s table method

1) Central differencing scheme O(`Delta x^2`)

`((del^2x)/(delx^2) =a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2))/(Delta x^2)`

The Taylor table can be obtained by expanding each function using the Taylor series

The Taylor Table for the central differencing scheme,

	f(x)	∆x*f\'(x)	(∆x^2)*f\'\'(x)	(∆x^3)*f\'\'\'(x)	(∆x^4)*f\'\'\'\'(x)	(∆x^5)*f\'\'\'\'\'(x)	(∆x^6)*f\'\'\'\'\'\'(x)
a*f(i-2)	a	(-2a)	4a/2	(-8a/6)	16a/24	(-32a/120)	64a/720
b*f(i-1)	b	(-b)	b/2	(-b/6)	b/24	(-b/120)	b/720
c*f(i)	c	0	0	0	0	0	0
d*f(i+1)	d	d	d/2	d/6	d/24	d/120	d/720
e*f(i+2)	e	2e	4e/2	8e/6	16e/24	32e/120	64e/720
∑	0	0	1	0	0	?	?

a,b,c,d,e can be obtained by converting the table into matrix form and then solve them using Matlab code.

2) Right side skewed scheme O(`Deltax^2`)

`(delx^2)/(delx^2) =(a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5))/(Deltax^2)`

The Taylor table can be obtained by expanding each function using the Taylor series

The Taylor table for right side skewed scheme,

	f(x)	∆x*f\'(x)	(∆x^2)*f\'\'(x)	(∆x^3)*f\'\'\'(x)	(∆x^4)*f\'\'\'\'(x)	(∆x^5)*f\'\'\'\'\'(x)	(∆x^6)*f\'\'\'\'\'\'(x)
a*f(i)	a	0	0	0	0	0	0
b*f(i+1)	b	b	b/2	b/6	b/24	b/120	b/720
c*f(i+2)	c	2c	4c/2	8c/6	16c/24	32c/120	64c/720
d*f(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120	729d/720
e*f(i+4)	e	4e	16e/2	64e/6	256e/24	1024e/120	4096e/720
g*f(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120	15625g/720
∑	0	0	1	0	0	0	?

a,b,c,d,e,g can be obtained by converting the table into matrix form and then solve them using Matlab code.

3) Left side skewed scheme O(`Deltax^2`)

`(del^2u)/(delx^2)=(a*f(i-5)+b*f(i-4)+c*f(i-3)+d*f(i-2)+e*f(i-1)+g*f(i))/(Deltax^2)`

The Taylor table can be obtained by expanding each function using the Taylor series

The Taylor table for left side skewed scheme,

	f(x)	∆x*f\'(x)	(∆x^2)*f\'\'(x)	(∆x^3)*f\'\'\'(x)	(∆x^4)*f\'\'\'\'(x)	(∆x^5)*f\'\'\'\'\'(x)	(∆x^6)*f\'\'\'\'\'\'(x)
a*f(i-5)	a	(-5a)	25a/2	(-125a/6)	625a/24	(-3125a/120)	15625a/720
b*f(i-4)	b	(-4b)	16b/2	(-64b/6)	256b/24	(-1024b/120)	4096b/720
c*f(i-3)	c	(-3c)	9c/2	(-27c/6)	81c/24	(-243c/120)	729c/720
d*f(i-2)	d	(-2d)	4d/2	(-8d/6)	16d/24	(-32d/120)	64d/720
e*f(i-1)	e	(-e)	e/2	(-e/6)	e/24	(-b/120)	b/720
g*f(i)	g	0	0	0	0	0	0
∑	0	0	1	0	0	0	?

a,b,c,d,e,g can be obtained by converting the table into matrix form and then solve them using Matlab code.

4)Forward differencing scheme O(`Deltax`)

`(del^2u)/(delx^2)=(f(i)-2*f(i+1)+f(i+2))/(Deltax^2)`

5)Backward differencing scheme O(`Deltax`)

`(del^2u)/(delx^2)=(f(i-2)-2*f(i+1)+f(i))/(Deltax^2)`

Matlab function for solving second order derivative using the center differencing scheme and also calculate the truncation error:

function out=fourth_order_center_approx(x,dx)
A=[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=inv(A)*B;
a=X(1);
b=X(2);
c=X(3);
d=X(4);
e=X(5);
analytical_value=-2*exp(x)*sin(x);
center_difference=((a*exc(x-(2*dx)))+(b*exc(x-dx))+(c*exc(x))+(d*exc(x+dx))+(e*exc(x+(2*dx))))/(dx^2);
out=abs(analytical_value-center_difference);
end

Matlab function for solving second order derivative using right side skewed scheme and also calculate the truncation error:

function out=right_side_skewed_approx(x,dx)
A=[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 4/3 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];
B=[0;0;1;0;0;0];
X=inv(A)*B;
a=X(1);
b=X(2);
c=X(3);
d=X(4);
e=X(5);
g=X(6);
analytical_value=-2*exp(x)*sin(x);
right_side_skewed_approx=((a*exc(x))+(b*exc(x+dx))+(c*exc(x+(2*dx)))+(d*exc(x+(3*dx)))+(e*exc(x+(4*dx)))+(g*exc(x+(5*dx))))/(dx^2);
out=abs(analytical_value-right_side_skewed_approx);
end

Matlab function for solving second order derivative using left side skewed scheme and also calculate the truncation error:

function out=left_side_skewed_approx(x,dx)
A=[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];
B=[0;0;1;0;0;0];
X=inv(A)*B;
a=X(1);
b=X(2);
c=X(3);
d=X(4);
e=X(5);
g=X(6);
analytical_value=-2*exp(x)*sin(x);
left_side_skewed_approx=((a*exc(x-(5*dx)))+(b*exc(x-(4*dx)))+(c*exc(x-(3*dx)))+(d*exc(x-(2*dx)))+(e*exc(x-dx))+(g*exc(x)))/(dx^2);
out=abs(analytical_value-left_side_skewed_approx);
end

Matlab function for solving second order derivative using forward differencing scheme(FDS) and also calculate the truncation error:

function out=first_order_fds_approx(x,dx)
analytical_value=-2*exp(x)*sin(x);
first_order_fds_approx=(exc(x)-(2*exc(x+dx))+exc(x+(2*dx)))/(dx^2);
out=abs(analytical_value-first_order_fds_approx);
end

Matlab function for solving second order derivative using backward differencing scheme(BDS) and also calculate the truncation error:

function out=first_order_bds_approx(x,dx)
analytical_value=-2*exp(x)*sin(x);
first_order_bds_approx=(exc(x)-(2*exc(x-dx))+exc(x-(2*dx)))/(dx^2);
out=abs(analytical_value-first_order_bds_approx);
end

For ease of coding, the separate function has been created for a given function f(x)=exp(x)*cos(x)

function a=exc(b)
a=exp(b)*cos(b);
end

Matlab code for plotting the errors vs dx

clear all
close all
clc
x=pi/3;
dx=linspace(pi/4,pi/40000,20);
for i=1:length(dx)
    center_diff_error(i)=fourth_order_center_approx(x,dx(i));
    right_side_error(i)=right_side_skewed_approx(x,dx(i));
    left_side_error(i)=left_side_skewed_approx(x,dx(i));
    fds_error(i)=first_order_fds_approx(x,dx(i));
    bds_error(i)=first_order_bds_approx(x,dx(i));
end
loglog(dx,center_diff_error,\'y\')
hold on
loglog(dx,right_side_error\',\'g\')
loglog(dx,left_side_error,\'b\')
loglog(dx,fds_error,\'k\')
loglog(dx,bds_error,\'r\')
legend(\'center differencing scheme\',\'right side skewed scheme\',\'left side skewed scheme\',\'FDS\',\'BDS\',\'location\',\'northwest\')
xlabel(\'dx\')
ylabel(\'Truncation error\')
title(\'Truncation error vs dx\')

Results:

1)Truncation errors vs dx

From the above graph, it was noticed that the left-sided scheme seemed to be unstable for a certain range of dx value. The central differencing scheme has a lower error value than other schemes.

2)The skewed scheme has a particular purpose where the central differencing scheme can\'t be performed. while solving the problem, in such a case the grid point doesn\'t have a neighbor grid points where skewed schemes can give an accurate solution.