Taylor table method and Matlab code - Computing higher order approximations using skewed and symmetric stencils

Objective:

1. To derive The following schemes for the 4th order approximations of the second order derivative
2. Central difference.
3. Skewed right sided difference.
4. Skewed left sided difference.
5. Forward difference scheme.
6. Backward difference scheme.

B. To provide a plot that compares the absolute error between the above mentioned schemes.

First let us choose the numerical stencils for the above mentioned schemes:

1. Central difference scheme:

In order to calculate the forth order approximation of a second order derivative five points are selected for a stencil. In CDS two points on each side of the function is selected. i.e., two on the left(-ve) and two on the right(+ve). The selected points are i-2 ,i-1, i, i+1, i+2. Therefore, The Taylor series for the CDS scheme can be written as:

In order to calculate constants for above scheme, Taylor’s table method is used:

	f(i)	∆xf’(i)	∆x2f’’(i)	∆ x3 f’’’(i)	∆x4 (i)	∆x5 (i)	∆x6 (i)
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	?	?

Writing down in matrix form:

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

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

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

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

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

The above equations are solved in matlab and the values for the constants a,b,c,d,e are found: a = -0.0833, b = 1.3333, c = -2.5000, d = 1.3333, e = -0.0833

1. Skewed right sided difference :

In order to calculate the forth order approximation of a second order derivative five points are selected for a stencil. In SRSD all points on right side of the function is selected. i.e., four points on the right(+ve) side of the function is selected.

The selected points are i, i+1, i+2,i+3,i+4.

Therefore, The Taylor series for the CDS scheme can be written as:

In order to calculate constants for above scheme, Taylor’s table method is used:

	f(i)	∆xf’(i)	∆x2f’’(i)	∆ x3 f’’’(i)	∆x4 (i)	∆x5 (i)	∆x6 (i)
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
	0	0	1	0	0	?	?

Writing down in matrix form:

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

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

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

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

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

The above equations are solved in matlab and the values for the constants a,b,c,d,e are found:

a = 1.3396, b = -2.9892, c = 1.9301, d = -0.2508, e = -0.0295

1. Skewed left sided difference :

In order to calculate the forth order approximation of a second order derivative five points are selected for a stencil. In SLSD all points on right side of the function is selected. i.e., four points on the right(-ve) side of the function is selected.

The selected points are i-4, i-3, i-2,i-1, i

Therefore, The Taylor series for the CDS scheme can be written as: In order to calculate constants for above scheme, Taylor’s table method is used:

	f(i)	∆xf’(i)	∆x2f’’(i)	∆ x3 f’’’(i)	∆x4 (i)	∆x5 (i)	∆x6 (i)
a*f(i-4)	a	-4a	(16a)/2	-(64a)/6	(256a)/24	-(1024a)/120	(4096a)/720
b*f(i-3)	b	-3b	(9b)/2	-(27b)/6	(81b)/24	-(243b)/120	(729b)/720
c*f(i-2)	C	-2c	(4c)/2	-(8c)/6	(16c)/24	-(32c)/120	(64c)/720
d*f(i-1)	d	-d	d/2	-d/6	d/24	-d/120	d/720
e*f(I)	e	0	0	0	0	0	0
	0	0	1	0	0	?	?

Writing down in matrix form:

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

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

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

-32a/3 - 9b/2 - 4c/3 - d/6 + 0 = 0 ..........(4)

32a/3 + 27b/8 + 2c/3 + d/24 + 0 = 0 ...........(5)

The above equations are solved in matlab and the values for the constants a,b,c,d,e are found:

a = 0.9166, b = -4.6666, c = 9.5000, d = -8.6666, e = 2.9166

The plot:

Plot 1(Comparison between CDS, SLDS, SRDS):

Plot 2(Comparison between CDS, SLDS, SRDS, FD and BD):

Following statements can be interpreted from above graphs:

1. The forward difference and backward difference scheme suits only for 1st order approx error of the first order derivative. Since it capture information only from one point (either right or left) it results in a lesser truncation error thus error value is high.

2. whereas in central or skewed scheme the information gathered is more because in order to calculate a fourth oder approximation of a second order derivative, The sycheme requires atleast five points around the function thus the error becomes low.

3. Among all the schemes central difference schemes hold the maximum accuracy with a lowest error value. This is because it gather informations from the point from both sides(symmetricity).

How is a skewed scheme useful?

Central difference scheme always requires equal no. of points on both side of the function which do not always go well, some systems wont be able to provide those symmetric points. Thus the skewed scheme always bemomes helpful when the numerical stencil is not symmetric and with skewed systems higher order accuracy can be obtained with given information without any constraints.

The matlab program is attached below:

The main program:

clear all
close all
clc

%calculating analytical derivative of function;
%given function, f(x)=exp(x).sin(x);
%f\'\'(x)=-2.exp(x).sin(x);

%%Assuming x=pi/3
x = pi/3;

%%Defining the grid size
dx=linspace(pi/10,pi/500,15);

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

%Calculating error
for i = 1:length(dx)
    error_central(i) =central_approximation(x,dx(i));
    error_right(i) = right_skewed_approximation(x,dx(i));
    error_left(i)=left_skewed_approximation(x,dx(i));
    error_forward(i)=forward_difference(x,dx(i));
    error_backward(i)=backward_difference(x,dx(i));
  
end

figure(1)
loglog(dx,error_central,\'color\', \'r\');
hold on
loglog(dx,error_right,\'color\', \'y\');
hold on
loglog(dx,error_left,\'color\', \'g\');

xlabel(\'dx\')
ylabel(\'truncation Error\')
title(\'Comparision between different schemes for a second order derivative of 4th order accuracy\')
legend(\'Central difference\',\'Skewed right sided difference\',\' Skewed left sided difference\')
%%Figure(2)
figure(2)
loglog(dx,error_central);
hold on
loglog(dx,error_right);
hold on
loglog(dx,error_left);
hold on
loglog(dx,error_forward);
hold on
loglog(dx,error_backward);
xlabel(\'dx\')
ylabel(\'truncation Error\')
title(\'second order derivative of 4th order accuracy vs FD vs BD\')
legend(\'CD\',\'SRSD\',\'SLSD\',\'Forward difference\',\'Backward difference\')

Program for function central_approximation:

function out = central_approximation(x,dx)

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

f = @(x) exp(x)*cos(x);

%calcuation of arbitary constants (a, b, c, d, e) for numerical derivative
syms a b c d e
eqn1 = a + b + c + d + e == 0;
eqn2 = -2*a - b + 0 + d + 2*e == 0;
eqn3 = 2*a + b/2 + 0 + d/2 + 2*e == 1;
eqn4 = -8*a/6 - b/6 + 0 + d/6 + 8*e/6 == 0;
eqn5 = 16*a/24 + b/24 + 0 + d/24 + 16*e/24 == 0;

%Solving Equations to calcualte constants
X = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [a, b, c, d, e]);

%Extracting constants
a = X.a;
b = X.b;
c = X.c;
d = X.d;
e = X.e;

%numerical approximation
central_differencing = (a*f(x-2*dx)+b*f(x-dx)+c*f(x)+d*f(x+dx)+e*f(x+2*dx))/(dx^2);

out = abs(central_differencing - analytical_derivative);

end

Program for function left_skewed_approximation:

function out = left_skewed_approximation(x,dx)

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

f = @(x) exp(x)*cos(x);

%calcuation of arbitary constants (a, b, c, d, e) for numerical derivative
syms a b c d e
eqn1 = a + b + c + d + e == 0;
eqn2 = -4*a - 3*b - 2*c - d == 0;
eqn3 = 16*a/2 + 9*b/2 + 4*c/2 + d/2 == 1;
eqn4 = -64*a/6 -27*b/6 - 8*c/6 -d/6 == 0;
eqn5 = 256*a/24 + 81*b/24 + 16*c/24 + d/24 == 0;

%Solving Equations to calcualte constants
X = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [a, b, c, d, e]);

%Extracting constants
a = X.a;
b = X.b;
c = X.c;
d = X.d;
e = X.e;

%numerical approximation
left_skewed_differencing = (a*f(x-4*dx)+b*f(x-3*dx)+c*f(x-2*dx)+d*f(x-dx)+e*f(x))/(dx^2);

out = abs(left_skewed_differencing - analytical_derivative);

end

Program for function right_skewed_approximation:

function out = right_skewed_approximation(x,dx)

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

f = @(x) exp(x)*cos(x);

%calcuation of arbitary constants (a, b, c, d, e) for numerical derivative
syms a b c d e
eqn1 = a + b + c + d + e == 0;
eqn2 = b + 2*c + 3*d + 4*e == 0;
eqn3 = b/2 + 4*c/2 + 9*d/2 + 16*e/2 == 1;
eqn4 = b/6 + 8*c/6 + 27*d/6 + 64*e/2 == 0;
eqn5 = b/24 + 16*c/24 + 81*d/24 + 256*e/24 == 0;

%Solving Equations to calcualte constants
X = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [a, b, c, d, e]);

%Extracting consytants
a = X.a;
b = X.b;
c = X.c;
d = X.d;
e = X.e;

%numerical approximation
right_skewed_differencing = (a*f(x)+b*f(x+dx)+c*f(x+2*dx)+d*f(x+3*dx)+e*f(x+4*dx))/(dx^2);

out = abs(right_skewed_differencing - analytical_derivative);

end

Function for forward difference scheme:

%%Function for forward difference
function out = forward_difference(x,dx)
%%Analytical value
analytical=-2*exp(x)*sin(x);
%%Forward difference
forward= (exp(x+dx)*cos(x+dx)-exp(x)*cos(x))/dx;
%Truncation Error
out = abs(forward-analytical);
end

Function for backward difference scheme:

%%%Function for backward difference
function out = backward_difference(x,dx)
%%Analytical value
analytical=-2*exp(x)*sin(x);
%%Forward difference
backward= (exp(x)*cos(x)-exp(x-dx)*cos(x-dx))/dx;
%Truncation Error
out = abs(backward-analytical);
end