Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

central difference scheme(CD):

no of nodes=p+q-1

where p=the order of accuracy

q=order of derivative

no of nodes for CD=4+2-1=5

Let the equation be

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

the taylor table is listed below:

`dx^2`f''(i)	f(i)	f'(i)dx	f''(i)`dx^2`	f'''(i) `dx^3`	f''''(i)`dx^4`
af(i-2)	a	-2a	2a	-4a/3	2a/3
b(fi-1)	b	-b	b/2	-b/6	b/24
cf(i)	c	0	0	0	0
df(i+1)	d	d	d/2	d/6	d/24
ef(i+2)	e	2e	2e	4e/3	2e/3
sum of co.eff	0	0	1	0	0

the matrix equation is given by,

`[[1,1,1,1,1],[-2,-1,0,1,2],[2,0.5,0,0.5,2],[-4/3,-1/6,0,1/6,4/3],[2/3,1/24,0,1/24,2/3]]**[[a],[b],[d],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

sollving this we get

a=-0.0833,b=1.3333,c=-2.5,d=1.3333,e=-0.0833

Right skewed scheme:

no of nodes=p+q=4+2=6

Let the equation be

`(del^2u)/(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))/(delx^2)`

the taylor table is listed below:

`dx^2`f''(i)	f(i)	f'(i)dx	f''(i)`dx^2`	f'''(i) `dx^3`	f''''(i)`dx^4`	f'''''(i)`dx^5`
af(i)	a	0	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24	b/120
cf(i+2)	c	2c	2c	4c/3	2c/3	4c/15
df(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120
ef(i+4)	e	4e	8e	32e/3	256e/24	1024e/120
gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120
sum of co.eff	0	0	1	0	0	0

solving the matrix equation we get,

`a=3.75,b=-12.833,c=17.833,d=-13,e=5.0833,g=-0.833`

Left skewed scheme:

no of nodes=p+q=4+2=6

Let the equation be

`(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))/(delx^2)`

the taylor table is listed below:

`dx^2`f''(i)	f(i)	f'(i)dx	f''(i)`dx^2`	f'''(i) `dx^3`	f''''(i)`dx^4`	f'''''(i)`dx^5`
af(i-5)	a	-5a	25a/2	-125a/6	625a/24	-3125a/120
bf(i-4)	b	-4b	8b	-32b/3	256b/24	-1024b/120
cf(i-3)	c	-3c	9c/2	-27c/6	81c/24	-243c/120
df(i-2)	d	-2d	2d	-4d/3	2d/3	-4d/15
ef(i-1)	e	-e	e/2	-e/6	e/24	-e/120
gf(i)	g	0	0	0	0	0
sum of co.eff	0	0	1	0	0	0

solving the matrix equation we get,

`a=-0.833,b=5.0833,c=-13,d=17.833,e=-12.833,g=3.75`

Matlab code:

clear all
close all
clc
x=pi/3;
n=20;
dx=linspace(pi/5,pi/500,n);
% analytical derivative
analytical_derivative=-2*exp(x)*sin(x);
for i=1:n
    central_diff(i)=cen_fourth_od(x,dx(i));
    right_skew(i)=right_sided(x,dx(i));
    left_skew(i)=left_sided(x,dx(i));
    error_cen(i)=abs(analytical_derivative-central_diff(i));
    error_rt(i)=abs(analytical_derivative-right_skew(i));
    error_lt(i)=abs(analytical_derivative-left_skew(i));
end

loglog(dx,error_cen,'r')
hold on
loglog(dx,error_rt,'b')
hold on
loglog(dx,error_lt,'m')
xlabel('range of dx')
ylabel('truncation error')
legend('central difference','right skewed','left skewed')

function [out]=cen_fourth_od(x,dx)
 [out]=((-0.0833)*exp(x-2*dx)*cos(x-2*dx)+(1.3333)*exp(x-dx)*cos(x-dx)-2.5*exp(x)*cos(x)+(1.3333)*exp(x+dx)*cos(x+dx)-(0.0833)*exp(x+2*dx)*cos(x+2*dx))/(dx^2);
end

function [out]=right_sided(x,dx)
    [out]=(3.75*exp(x)*cos(x)-12.8333*exp(x+dx)*cos(x+dx)+17.8333*exp(x+2*dx)*cos(x+2*dx)-13*exp(x+3*dx)*cos(x+3*dx)+5.0833*exp(x+4*dx)*cos(x+4*dx)-0.8333*exp(x+5*dx)*cos(x+5*dx))/(dx^2);
end

function [out]=left_sided(x,dx)
    [out]=((-0.8333*exp(x-5*dx)*cos(x-5*dx)+5.0833*exp(x-4*dx)*cos(x-4*dx))-13*exp(x-3*dx)*cos(x-3*dx)+17.8333*exp(x-2*dx)*cos(x-2*dx)-12.8333*exp(x-dx)*cos(x-dx)+3.75*exp(x)*cos(x))/(dx^2);
end

Results and discussions:

as we see from the graph that the error in central differencing scheme is less compared to one sided schemes.

for computing central difference information from both the sides of point to be considered is necessary. whereas in case of one sided scheme information from either one side is enough

and thus CD cannot be employed for computing values at end nodes but one sided scheme can be employed