Taylor table method and Matlab code

Objective:

To derive the fourth order approximations of the second order derivative using central differencing, skewed right difference and skewed left difference.

1) Central differencing.

∂²f/∂x² ≈af(i-2) + bf(i-1) +cf(i) + df(i+1) + ef(i+2) (1)

Using the Taylor\'s Table method to find out the co-efficient of a,b,c,d,e.

	f(i)	f’(i)* Δx	f’’(i)* Δx2	f’’’(i)* Δx3	fiv(i)* Δx4	fv(i)* Δx5	fvi(i)* Δx6
a*f(i-2)	a	-2a	(4/2)a	(-8/6)a	(16/24)a	(-32/120)a	(64/240)a
b*f(i-1)	b	-b	(1/2)b	(-1/6)b	(1/24)b	(-1/120)b	(1/720)b
c*f(i)	c	0	0	0	0	0	0
d*f(i+1)	d	d	(1/2)d	(1/6)d	(1/24)d	(1/120)d	(1/240)d
e*f(i+2)	e	2e	(4/2)e	(8/6)e	(16/24)e	(32/120)e	(64/240)e

There are five unknowns and hence five equations are needed. From the above table we can create five equations which are as follows:

a+b+c+d+e=0

-2a-b+0+d+2e=0

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

(-8/6)a-(1/6)b+0+(1/6)d+(8/6)e=0

(16/24)a+(1/24)b+0+(1/24)d+(16/24)e=0

On solving these equations we get the values of a,b,c,d and e as

a = -0.1786

b = 1.7143

c = -3.0714

d = 1.7143

e = -0.1786

Now the equation tends to

∂2f/∂x2 = a f(i-2)+b f(i-1)+c f(i)+d f(i+1)+e f(i+2) =

(2a+0.5b+0+0.5d+2e) f’’(i) Δx² +((64/240)a+(1/240)b+0+(1/240)d+(64/240)e) f^vi(i) Δx⁶

Dividing the above equation by Δx² , we get

[-0.1786*f(i-2)+1.7143*f(i-1)-3.0714*f(i)+1.7143*f(i+1)-0.1786*f(i+2)]/ Δx² =

(2a+0.5b+0+0.5d+2e) f’’(i) +((64/240)a+(1/240)b+0+(1/240)d+(64/240)e) f^vi(i) Δx⁴

Now the truncation error term denotes that the order of the approximation is 4th order.

2) Skewed right difference.

∂²f/∂x² ≈af(i) + bf(i+1) +cf(i+2) + df(i+3) + ef(i+4) (2)

Using the Taylor\'s Table method to find out the co-efficient of a,b,c,d,e.

	f(i)	f’(i)* Δx	f’’(i)* Δx2	f’’’(i)* Δx3	fiv(i)* Δx4	fv(i)* Δx5	fvi(i)* Δx6
a*f(i)	a	0	0	0	0	0	0
b*f(i+1)	b	b	(1/2)b	(1/6)b	(1/24)b	(1/120)b	(1/720)b
c*f(i+2)	c	2c	2c	(8/6)c	(16/24)c	(32/120)c	(64/720)c
d*f(i+3)	d	3d	(9/2)d	(27/6)d	(81/24)d	(243/120)d	(729/720)d
e*f(i+4)	e	4e	(16/2)e	(64/6)e	(256/24)e	(1024/120)e	(4096/720)e

There are five unknowns and hence five equations are needed. From the above table we can create five equations which are as follows:

a+b+c+d+e=0

0+b+2c+3d+4e=0

0+(1/2)b+2c+(9/2)d+(16/2)e=1

0+(1/6)b+(8/6)c+(27/6)d+(64/6)e=0

0+(1/24)b+(16/24)c+(81/24)d+(256/24)e=0

On solving these equations we get the values of a,b,c,d and e as

a = 2.9167

b = -8.6667

c = 9.5

d = -4.6667

e = 0.9167

Now the equation tends to

∂2f/∂x2 = a f(i)+b f(i+1)+c f(i+2)+d f(i+3)+e f(i+4) =

(0+0.5b+2c+(9/2)d+(16/2)e) f’’(i) Δx² +(0 +(1/720)b+(64/720)c+(729/720)d+(4096/720)e) f^vi(i) Δx⁶

Dividing the above equation by Δx² , we get

[2.9167*f(i)-8.6667*f(i+1)+9.5*f(i+2)-4.6667*f(i+3)+0.9167*f(i+4)]/ Δx² =

(0+0.5b+2c+(9/2)d+(16/2)e) f’’(i) + (0+(1/720)b+(64/720)c+(729/720)d+(4096/720)e) f^vi(i) Δx⁴

Now the truncation error term denotes that the order of the approximation is 4th order.

3) Skewed left differencing.

∂²f/∂x² ≈af(i-4) + bf(i-3) +cf(i-2) + df(i-1) + ef(i) (3)

Using the Taylor\'s Table method to find out the co-efficient of a,b,c,d,e.

	f(i)	f’(i)* Δx	f’’(i)* Δx2	f’’’(i)* Δx3	fiv(i)* Δx4	fv(i)* Δx5	fvi(i)* Δx6
a*f(i-4)	a	-4a	8a	-(64/6)a	(256/24)a	-(1024/120)a	(4096/720)a
b*f(i-3)	b	-3b	(9/2)b	-(27/6)b	(81/24)b	-(243/120)b	(729/720)b
c*f(i-2)	c	-2c	2c	-(8/6)c	(16/24)c	-(32/120)c	(64/720)c
d*f(i-1)	d	-d	(1/2)d	-(1/6)d	(1/24)d	-(1/120)d	(1/720)d
e*f(i)	e	0	0	0	0	0	0

There are five unknowns and hence five equations are needed. From the above table we can create five equations which are as follows:

a+b+c+d+e=0

-4a-3b-2c-d+0=0

8a+(9/2)b+2c+(1/2)d+0=1

-(64/6)a-(27/6)b-(8/6)c-(1/6)d+0=0

(256/24)a+(81/24)b+(16/24)c+(1/24)d+0=0

On solving these equations we get the values of a,b,c,d and e as

a = 0.9167

b = -4.6667

c = 9.5

d = -8.6667

e = 2.9167

Now the equation tends to

∂2f/∂x2 = a f(i)+b f(i+1)+c f(i+2)+d f(i+3)+e f(i+4) =

(8a+4.5b+2c+0.5d+0) f’’(i) Δx² +((4096/720)a+(729/720)b+(64/720)c+(1/720)d+0) f^vi(i) Δx⁶

Dividing the above equation by Δx² , we get

[0.9167*f(i-4)-4.6667*f(i-3)+9.5*f(i-2)-8.6667*f(i-1)+2.9167*f(i)]/ Δx² =

(8a+4.5b+2c+0.5d+0) f’’(i)

+((4096/720)a+(729/720)b+(64/720)c+(1/720)d+0) f^vi(i) Δx⁴

Now the truncation error term denotes that the order of the approximation is 4th order.

% second order derivative of the given function

clear all
close all
clc

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

x=pi/3;
dx=pi/40;

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

% numerical derivative of the given function
% fourth order approximation using central differencing

central_differencing=(-2.1432*exp(x-2*dx)*cos(x-2*dx)+20.5716*exp(x-dx)*cos(x-dx)-36.8568*exp(x)*cos(x)+20.5716*exp(x+dx)*cos(x+dx)-2.1432*exp(x+2*dx)*cos(x+dx))/(12*dx^2);

% skewed right side differencing
skewed_right_side_differencing= (35*exp(x)*cos(x)-104*exp(x+dx)*cos(x+dx)+114*exp(x+2*dx)*cos(x+2*dx)-56*exp(x+3*dx)*cos(x+3*dx)+11*exp(x+4*dx)*cos(x+4*dx))/(12*dx^2);

% skewed left side diffrencing
skewed_left_side_differencing=(11*exp(x-4*dx)*cos(x-4*dx)-56*exp(x-3*dx)*cos(x-3*dx)+114*exp(x-2*dx)*cos(x-2*dx)-104*exp(x-dx)*cos(x-dx)+35*exp(x)*cos(x))/(12*dx^2);

% error

error_central_differencing=abs(central_differencing-analytical_derivative);
error_skewed_right_side_differencing=abs(skewed_right_side_differencing-analytical_derivative);
error_skewed_left_side_differencing=abs(skewed_left_side_differencing-analytical_derivative);

%comparison with FDS and BDS Methods

forward_difference=((exp(x+dx)*cos(x+dx))-(exp(x)*cos(x)))/dx;
backward_difference=((exp(x-dx)*cos(x-dx))-(exp(x)*cos(x)))/dx;

error_fds= abs(forward_difference-analytical_derivative)
error_bds=abs(backward_difference-analytical_derivative)

fname={\'central differencing\',\'skewed right \',\'skewed left\',\'forward difference\',\'backward difference\'}
figure(1)
x=[1:5];
y=[error_central_differencing error_skewed_right_side_differencing error_skewed_left_side_differencing error_fds error_bds];

bar(x,y);
set(gca,\'xTick\',1:length(fname),\'xTickLabel\',fname);
xlabel(\'scheme\')
ylabel(\'truncation error\')

From the above code, the absolute error between analytical derivative and the 3 schemes have been developed. It also shows the comparison between Forward and Backward Differencing schemes.

function out=central_differencing(x,dx)


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

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

% numerical derivative
central_differencing=(-2.1432*exp(x-2*dx)*cos(x-2*dx)+20.5716*exp(x-dx)*cos(x-dx)-36.8568*exp(x)*cos(x)+20.5716*exp(x+dx)*cos(x+dx)-2.1432*exp(x+2*dx)*cos(x+dx))/(12*dx^2);

out=abs(central_differencing-analytical_derivative);
end

function out=skewed_right_side_differencing(x,dx)


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

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

% numerical derivative
skewed_right_side_differencing=(35*exp(x)*cos(x)-104*exp(x+dx)*cos(x+dx)+114*exp(x+2*dx)*cos(x+2*dx)-56*exp(x+3*dx)*cos(x+3*dx)+11*exp(x+4*dx)*cos(x+4*dx))/(12*dx^2);

out=abs(skewed_right_side_differencing-analytical_derivative);
end

function out=skewed_left_side_differencing(x,dx)


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

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

% numerical derivative
skewed_left_side_differencing=(11*exp(x-4*dx)*cos(x-4*dx)-56*exp(x-3*dx)*cos(x-3*dx)+114*exp(x-2*dx)*cos(x-2*dx)-104*exp(x-dx)*cos(x-dx)+35*exp(x)*cos(x))/(12*dx^2);

out=abs(skewed_left_side_differencing-analytical_derivative);
end

clear all
close all
clc
x=pi/3;
dx=linspace(pi/4,pi/400,5);
for i= 1: length(dx)

error_central_difference(i)= central_differencing(x,dx(i));
error_skewed_right_sided_differencing(i)= skewed_right_side_differencing(x,dx(i));
error_skewed_left_sided_differencing(i)=skewed_left_side_differencing(x,dx(i));
end 

figure(1)
loglog(dx,error_central_difference)
hold on
loglog(dx,error_skewed_right_sided_differencing,\'r\')
hold on
loglog(dx,error_skewed_left_sided_differencing,\'g\')

xlabel(\'step size\')
ylabel(\'truncation error\')
legend(\'CD \',\'Skewed right\',\'skewed left\')

The above bar graph shows the comparison between the different schemes and the error.

Since the skewed errors are very small compared to the central differencing scheme, it has not been highlighted.

The above log plot compares the absolute error with the above-mentioned schemes.

From the above plot, it is clear that as the step size increases, the error also increases. So a smaller step size is recommended.

The skewed schemes can be used anywhere on the set of available discrete data.

Suppose if we have discontinuities or if there is an imbalance in the number of points available on the grid, then we can use skewed schemes but not central differencing schemes.