Taylor table method and matlab code for a fourth order approximation of second order derivative

Aim: A)To derive the 4th order approximations of the second order derivative

1. Central difference
2. Skewed right sided difference
3. skewed left handed difference

B) To write a Matlab Program to evaluate the second order derivative of analytical function exp(x)*cos(x) and compare it with numerical approximations, and also to compare the absolute errors between the second order schemes.

C) second order schemes comparison with FDS and BDS.

1. Fourth order Central difference of the second derivative:

   `frac{∂^2 f}{∂x^2 }` = 4 th order accurate approximation using i -2, i -1, i, i +1, i +2

   Here, the stencil is symmetric, it is a symmetric 2nd order scheme.

   Using Taylor series method,

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

   Writing the tabular form, we have

	f(i)	Δx fI(i)	Δx2 fII(i)	Δx3 fIII(i)	Δx4 fIV(i)	Δx5 fV(i)	Δx6 fVI(i)
af(i-2)	a	-2a	(4a)/2	-(8a)/6	(16a)/24	-(32a)/120	(64a)/720
bf(i-1)	b	-b	b/2	-b/6	b/24	-b/120	b/720
cfi	c	0	0	0	0	0	0
df(i+1)	d	d	d/2	d/6	d/24	d/120	d/720
ef(i+2)	e	2e	(4e)/2	(8e)/6	(16e)/24	(32e)/120	(64e)/720
	0	0	1	0	0	?	?

On further simplification, we get

`[[1,1,1,1,1],[-2,-1,0,1,2],[4/2,1/2,0,1/2,4/2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]][(a),(b),(c),(d),(e),(f)] = [(0),(0),(1),(0),(0)]`

This is in the form of \'\'AX = B\'\' . Rewriting them in equation form,

a + b + c + d + e =0 ...........(2)

-2a - b + 0 + d + 2e = 0 ...........(3)

2a + b/2 + 0 + d/2 + 2e = 1 ...........(4)

-4/3 - b/6 + 0 + d/6 + 4/3 =0 ...........(5)

2a/3 + b/24 + 0 + d/24 + 2a/3 = 0 .....(6)

we have 5 equations and 5 unknowns, finally we get

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

substituting a,b,c,d,e in equation (i) so we get

`(∂^2 f)/(∂x^2 ) approx (-0.0833f(i-2) + 1.3333f(i-1) - 2.5f(i) + 1.3333f(i+1)+ 0.08333f(i+2))/(Δx^2)`

2. Fourth order skewed right sided difference of the second derivative:

`frac{∂^2 f}{∂x^2 }` = 4 th order accurate approximation using i , i +1, i, i +2, i +3,

i +4, i+5

Here, the stencil is not symmetric, but it is a skewed right sided difference of 2nd order scheme.

Using Taylor series method,

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

Writing the tabular form, we have

	f(i)	Δx fI(i)	Δx2 fII(i)	Δx3 fIII(i)	Δx4 fIV(i)	Δx5 fV(i)	Δx6 fVI(i)
af(i)	a	0	0	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24	b/120	b/720
cf(i+2)	c	2c	2c	(8c)/6	(16c)/24	(32c)/120	(64c)/720
df(i+3)	d	3d	(9d)/2	(27d)/6	(81d)/24	(243d)/120	(729d)/720
ef(i+4)	e	4e	(16e)/2	(32e)/3	(256e)/24	(1024e)/120	(4096e)/720
gf(i+5)	g	5g	(25g)/2	(125g)/6	(625g)/24	(3125g)/120	(15625g)/720
	0	0	1	0	0	?	?

On further simplification, we get

`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,2,9/2,8,25/2],[0,1/6,4/3,9/2,32/3,125/6],[0,1/24,2/3,27/8,32/3,625/24],[0,1/120,32/120,243/120,1024/120,3125/120]][(a),(b),(c),(d),(e),(g)] = [(0),(0),(1),(0),(0),(0)]`

This is in the form of \'\'AX = B\'\' . Rewriting them in equation form,

a + b + c + d + e +g=0 ...........(7)

0 + b + 2c + 3d + 4e + 5g= 0 ...........(8)

0 + b/2 + 2c + 9d/2 + 8e + 25g/2 = 1 ...........(9)

0 + b/6 + 4c/3 + 9d/2 + 32e/3 + 125g/6 =0 ...........(10)

0 + b/24 + 2c/3 + 27d/8 + 32e/3 + 625g/24 = 0 .....(11)

0 + b/120 + 32c/120+ 243d/120 + 1024e/120 + 3125g/120= 0.......(12)

we have 6 equations and 6 unknowns, finally we get

`a= 3.75, b= -12.8333, c= 17.8333, d= -13, e = -5.0833,g=-0.8333`

substituting a,b,c,d,e in equation (i) so we get

`(∂^2 f)/(∂x^2 ) approx (1.3396f(i) - 2.9892f(i+1) + 1.9301f(i+2) - 0.2508f(i+3) - 0.0295f(i+4))/(Δx^2)`

3. Fourth order skewed left sided difference of the second derivative:

`frac{∂^2 f}{∂x^2}` = 4 th order accurate approximation using i, i -1, i-2, i -3, i-4 and

i-5

Here, the stencil is not symmetric, but it is a skewed left sided difference of 2nd order scheme.

Using Taylor series method,

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

Writing the tabular form, we have

	f(i)	Δx fI(i)	Δx2 fII(i)	Δx3 fIII(i)	Δx4 fIV(i)	Δx5 fV(i)	Δx6 fVI(i)
af(i)	a	0	0	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24	b/120	b/720
cf(i+2)	c	2c	2c	(8c)/6	(16c)/24	(32c)/120	(64c)/720
df(i+3)	d	3d	(9d)/2	(27d)/6	(81d)/24	(243d)/120	(729d)/720
ef(i+4)	e	4e	(16e)/2	(32e)/3	(256e)/24	(1024e)/120	(4096e)/720
gf(i+5)	g	5g	(25g)/2	(125g)/6	(625g)/24	(3125g)/120	(15625g)/720
	0	0	1	0	0	?	?

On further simplification, we get

`[[1,1,1,1,1,1],[0,-1,-2,-3,-4,-5],[0,1/2,2,9/2,8,25/2],[0,-1/6,-4/3,-9/2,-32/3,-125/6],[0,1/24,2/3,27/8,32/3,625/24],[0,-1/120,-32/120,-243/120,-1024/120,-3125/120]][(a),(b),(c),(d),(e),(g)] = [(0),(0),(1),(0),(0),(0)]`

This is in the form of \'\'AX = B\'\' . Rewriting them in equation form,

`a + b + c + d + e + g=0 ...........(13)`

`0- b - 2c - 3d - 4e - 5g = 0 ..........(14)`

`0+ b/2 + 2c + {9d}/2 + 8e + {25g}/2= 1 ...........(15)`

`0 - b/6 - {4c}/3 - {9d}/2 - {32e}/3 - {125g}/6 =0 ...........(16)`

`0+ b/24 + {2c}/3 + {27d}/8 + {32e}/3 + {625g}/24 = 0 .....(17)`

`0 - b/120 - {32c}/120- {243d}/120 - {1024e}/120 - {3125g}/120 = 0.....(18))`

we have 5 equations and 5 unknowns, finally we get

`a= 3.75, b= -12.8333, c= 17.8333, d= -13, e = -5.0833,g=-0.8333`

Calling a function of fourth order second derivative Central differencing :

function out =fourth_order_second_derivative_central_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

exact_second_derivative = -2*exp(x)*sin(x);

% fourth order second derivative central differencing = (a*f(i-2) + b*f(i-1) + c*f(i) + d*f(i+1) + e*f(i+2))/(dx^2);
% The constants calculated using 5 equations are:  a = -0.0833, b = 1.3333, c = -2.5, d = 1.3333, e = -0.0833
fourth_order_second_derivative_central_differencing = ((-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)

% fourth order second derivation approximation(central differencing)
out = abs(fourth_order_second_derivative_central_differencing-exact_second_derivative)

end

Calling a function of fourth order second derivative Right side Skeweddifferencing:

function out =fourth_order_second_derivative_skewed_right_side_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

exact_second_derivative = -2*exp(x)*sin(x);

% numerical derivative
% skewed_right_sided_differencing = (a*f(i) + b*f(i+1) + c*f(i+2) + d*f(i+3) + e*f(i+4) + g*f(i+5)/(dx^2);
% The constants calculated using 5 equations are:  a = 3.75, b = -12.8333, c = 17.8333, d = -13.0000, e = 5.0833, g = -0.8333 
fourth_order_second_derivative_skewed_right_side_differencing = ((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.83333*exp(x+(5*dx))*cos(x+(5*dx))))/(dx^2)

% fourth order second derivation approximation(right sided differencing)
out = abs(fourth_order_second_derivative_skewed_right_side_differencing-exact_second_derivative)

end

Calling a function of fourth order second derivative Left side Skewed differencing:

function out =fourth_order_second_derivative_skewed_left_side_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

exact_second_derivative = -2*exp(x)*sin(x);

% numerical derivative
% skewed_left_sided_differencing = (a*f(i) + b*f(i-1) + c*f(i-2) + d*f(i-3) + e*f(i-4) + gf(i-5)/(dx^2);
% The constants calculated using 5 equations are:  a = 3.75, b = -12.8333, c = 17.8333, d = -13, e = 5.0833, g = -0.8333 
fourth_order_second_derivative_skewed_left_side_differencing = ((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)

% fourth order second derivation approximation(left sided differencing)
out = abs(fourth_order_second_derivative_skewed_left_side_differencing-exact_second_derivative)

end

Calling a function of forward Differencing:

function out = forward_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

exact_second_derivative = -2*exp(x)*sin(x);

% numerical derivative
% forward differencing = (f(x+dx)-f(x))/dx

forward_differencing = ((exp(x+dx)*cos(x+dx)) - (exp(x)*cos(x)))/dx

% First order approximation
out = abs(forward_differencing - exact_second_derivative)

end

Calling a function of Backward Differencing:

function out = backward_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

exact_second_derivative = -2*exp(x)*sin(x);

% numerical derivative
% forward differencing = (f(x+dx)-f(x))/dx

backward_differencing = ((exp(x)*cos(x)) - (exp(x-dx)*cos(x-dx)))/dx

% First order approximation
out = abs(backward_differencing - exact_second_derivative)

end

The Main Matlab Program:

clear all
close all 
clc

% analytical function = exp(x)*cos(x)
% f\"(x) = -2*exp(x)*sin(x)

% Assuming the input value of x
x = pi/3;
% Defining the grid size,dx
dx = linspace(pi/10,pi/100,10)

% second_order_derivative = -2*exp(x)*sin(x);
for i = 1:length(dx)

fourth_order_second_derivative_central_differencing_error(i) = fourth_order_second_derivative_central_differencing(x,dx(i))
fourth_order_second_derivative_skewed_right_side_differencing_error(i) = fourth_order_second_derivative_skewed_right_side_differencing(x,dx(i))
fourth_order_second_derivative_skewed_left_side_differencing_error(i) = fourth_order_second_derivative_skewed_left_side_differencing(x,dx(i))
forward_differencing_error(i) = forward_differencing(x,dx(i))
backward_differencing_error(i) = backward_differencing(x,dx(i))
end

figure(1)
loglog(dx,fourth_order_second_derivative_central_differencing_error)
hold on
loglog(dx, fourth_order_second_derivative_skewed_right_side_differencing_error,\'r\')
hold on
loglog(dx, fourth_order_second_derivative_skewed_left_side_differencing_error,\'g\')
xlabel(\'Grid size (dx)\')
ylabel(\'Error\')
legend(\'central differencing\',\'Skewed Right sided differencing\',\'skewed Left sided differencing\')
title(\'Comparison of errors for different fourth order second derivative schemes\')

figure(2)
loglog(dx,fourth_order_second_derivative_central_differencing_error)
hold on
loglog(dx, fourth_order_second_derivative_skewed_right_side_differencing_error,\'r\')
hold on
loglog(dx, fourth_order_second_derivative_skewed_left_side_differencing_error,\'g\')
hold on
loglog(dx, forward_differencing_error,\'y\')
hold on
loglog(dx, backward_differencing_error,\'black\')
xlabel(\'Grid size (dx)\')
ylabel(\'Error\')
legend(\'central differencing\',\'Skewed Right sided differencing\',\'skewed Left sided differencing\',\'forward differencing\',\'backward differencing\')
title(\'Comparison of errors for different fourth order second derivative schemes VS Forward Differencing VS Backward Differencing\')

Plot 1:

Plot 2:

The absolute error for skewed schemes and cenral scheme with many neighbouring points results is approximately a curve with respect to grid size dx. This is due to reason that the influence of many neighbouring points and a very small value of dx, generates large no.of airthmetic operations, which results in round off error. where for the same range of values dx, the error obtained by FDS and BDS are not curved, due to less no. of airthmetic operations.

Central differencing scheme error is less, but skewed schemes are mostly used because it can compute a single or higher order derivative at any corner node, in this case is no more central scheme is useful. If scheme is not symmetric in those cases only skewed schemes are helpful, which means no. of node points on right and left are not same .