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

Aim: To evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations obtained using skewed right side, skewed left side and central differencing fourth order approximations of second derivative.

Concept:

1)Skewed right sided difference:

Formula:`(del^2y)/(delx^2)=(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/h^2`

No. of nodes=P+Q=4+2=6

Taylor's table

	f(x)	hf'(x)	`h^2f''(x)`	`h^3f'''(x)`	`h^4f''''(x)`	`h^5f'''''(x)`
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	4c/2	8c/6	16c/24	32c/120
df(i+4)	d	3d	9d/2	27d/6	81d/24	243d/120
ef(i+4)	e	4e	16e/2	64e/6	256e/24	1024e/120
gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120
af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4) gf(i+5)	0	0	1	0	0	0

And then the 6 equations with 6 unknowns were obtained, expressed as AX=B. Thus X= inv(A)*B

a=3.7500
b=-12.8333
c= 17.8333
d= -13.0000
e= 5.0833
g= -0.8333

2)Skewed left sided difference:

Formula:`(del^2y)/(delx^2)=(af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i))/h^2`

No. of nodes=P+Q=4+2=6

Taylor's table

	f(x)	hf'(x)	`h^2f''(x)`	`h^3f'''(x)`	`h^4f''''(x)`	`h^5f'''''(x)`
af(i-5)	a	-5a	25a/2	-125a/6	625a/24	-3125a/120
bf(i-4)	b	-4b	16b/2	-64b/6	256b/24	-1024b/120
cf(i-3)	c	-3c	9c/2	-27c/6	81c/24	-243c/120
df(i-2)	d	-2d	4d/2	-8d/6	16d/24	-32d/120
ef(i-1)	e	-e	e/2	-e/6	e/24	-e/120
gf(i)	g	0	0	0	0	0
af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)	0	0	1	0	0	0

And then the 6 equations with 6 unknowns were obtained, expressed as CY=D Thus Y= inv(C)*D

a=-0.8333
b=5.0833
c=-13.0000
d=17.8333
e=-12.8333
g=3.7500

3)Central Difference:

Formula:`(del^2y)/(delx^2)=(af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/h^2`

No. of nodes=P+Q-1=4+2-1=5`

Taylor's table

	f(x)	hf'(x)	`h^2f''(x)`	`h^3f'''(x)`	`h^4f''''(x)`
af(i-2)	a	-2a	4a/2	-8a/6	16a/24
bf(i-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	4e/2	8e/6	16e/24
af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)	0	0	1	0	0

And then the 5 equations with 5 unknowns were obtained, expressed as EZ=F Thus Z= inv(E)*F.

a=-0.0833
b=1.3333
c=-2.5000
d=1.3333
e=-0.0833

In left and right sided difference, the coefficient of the term `h^6f''''''(x)`is non zero. This upon division by `h^2`leaves with an error in the fourth order.

Similar is the case with central difference too. The coefficients of `h^5f'''''(x)`must be zero even though that constrained was not directly placed while solving 5 equations in EZ=F.

All the coefficients were obtained by proper codes in MATLAB and were used to calculate the values of derivatives

Main Program code:

% PROGRAM FOR FINDING THE FOURTH ORDER APPROXIMATION OF SECOND DERIVATIVE
% OF EXP(X)*COSX NUMERICALLY AND ANALYTICALLY AND COMPARING ERRORS
clear all
close all
clc
%% DEFINING INPUTS AND ANALYTICAL DERIVATIVE OF EXP(X)*COS(X)
x=pi/3;
dx=linspace(pi/20,pi/400,20);
d2y_dx2=-2*exp(x)*sin(x);

%% OBTAINING SKEWED RIGHT SIDED DIFFERENCE COEFFICIENTS X FROM AX=B USING TAYLOR'S TABLE METHOD
A=[1 1 1 1 1 1;
   0 1 2 3 4 5;
   0 0.5 2 4.5 8 25/2
   0 1/6 8/6 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

%% OBTAINING SKEWED LEFT SIDED DIFFERENCE COEFFICIENTS Y FROM CY=D USING TAYLOR'S TABLE METHOD
C=[1 1 1 1 1 1;
   -5 -4 -3 -2 -1 0;
   25/2 16/2 9/2 4/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];
D=[0;0;1;0;0;0];
Y=inv(C)*D
%% OBTAINING CENTRAL DIFFERENCE COEFFICIENTS Z FROM EZ=F USING TAYLOR'S TABLE METHOD
E=[1 1 1 1 1;
   -2 -1 0 1 2;
   2 1/2 0 1/2 2;
   -4/3 -1/6 0 1/6 4/3;
   2/3 1/24 0 1/24 2/3];
F=[0;0;1;0;0];
Z=inv(E)*F
%% COMPUTING RIGHT SIDED SKEWED APPROXIMATED VALUE AND ERROR FROM TRUE VALUE
right_skd_value=(1./dx.^2).*(X(1).*exp(x).*cos(x)+X(2).*exp(x+dx).*cos(x+dx)+X(3).*exp(x+2*dx).*cos(x+2*dx)+X(4).*exp(x+3*dx).*cos(x+3*dx)+X(5).*exp(x+4*dx).*cos(x+4*dx)+X(6).*exp(x+5*dx).*cos(x+5*dx))
right_skd_error=abs(right_skd_value-d2y_dx2);
%% COMPUTING LEFT SIDED SKEWED APPROXIMATED VALUE AND ERROR FROM TRUE VALUE
left_skd_value=(1./dx.^2).*(Y(1).*exp(x-5*dx).*cos(x-5*dx)+Y(2).*exp(x-4*dx).*cos(x-4*dx)+Y(3).*exp(x-3*dx).*cos(x-3*dx)+Y(4).*exp(x-2*dx).*cos(x-2*dx)+Y(5).*exp(x-dx).*cos(x-dx)+Y(6).*exp(x).*cos(x))
left_skd_error=abs(left_skd_value-d2y_dx2);
%% COMPUTING CENTRAL SKEWED APPROXIMATED VALUE AND ERROR FROM TRUE VALUE
central_skd_value=(1./dx.^2).*(Z(1).*exp(x-2*dx).*cos(x-2*dx)+Z(2).*exp(x-dx).*cos(x-dx)+Z(3).*exp(x).*cos(x)+Z(4).*exp(x+dx).*cos(x+dx)+Z(5).*exp(x+2*dx).*cos(x+2*dx))
central_skd_error=abs(central_skd_value-d2y_dx2);
%% PLOTTING LOG LOG GRAPH FOR ERROR VS SPACING FOR RIGHT SIDED SKEWED, LEFT SIDED SKEWED AND CENTRAL SKEWED SECOND DIFFERENCES
loglog(dx,right_skd_error,'r');
hold on
loglog(dx,left_skd_error,'g');
loglog(dx,central_skd_error,'m')
legend('right skewed','left skewed','central skewed')
xlabel('ln(dx)')
ylabel('ln(error)')
title('ln(spacing) vs ln(error)')

1)The values of analytical derivative of the given function value was calculated at x=pi/3 and stored in the variable d2y_dx2

2)The value for spacing dx is chosen arbitrarily

3)The value of the coefficients in the expressions were calculated using the Taylor's table and solving the equations using matrix method

4)Using those coefficients obtained, the value of skewed differences-left and right and central differences were computed for different values of dx resulting in an array of values being obtained having same size as that of dx array .

5) Absolute values of error computed are clearly in the order of (spacing)^4

6)Then, the graph was plotted on a log-log scale to draw comparison between different numerical schemes.

Inference:

1)The central differencing scheme was highly accurate compared to other two schemes. Skewed left or skewed right schemes were more or less same with skewed left being slightly more accurate

2)But, the problem with central differencing scheme is that it takes information from both the directions to compute the value . So, there might be problem when the values are computed near the boundary. So, an alternative is to use skewed schemes to use only near boundaries and use central differencing schemes at all other points. But, the error incurred due to thiis going to affect all the values

3)Another alternative is to use Skewed differencing scheme everywhere. This is more common as it keeps things simpler without unnecessary complications.