Week 3 - Taylor table method and Matlab code

A REPORT ON TAYLOR TABLE

By the taylor series, it became possible to convert derivative functions into numerical algebraic functions.

And by comparing Numerical algebraic equation of derivative function(1^st order,2^nd order….4^th order….) by its analytical derivative, we get the values of error which is called as truncation error.

The error value contain grid spacing term(h or dx) and the power of h decides the accuracy.

Higher the power, higher will be the accuracy, hence lower the error.

Now,

f’(x) = f(x+h) –f(x)/(h) + O(h)

This first derivative by forward differencing is 1st order accurate having two nodes.

Similarly

f’(x) = f(x+h) – f(x-h)/(2h) + O(h^2)

This first derivative by central differencing is 2^nd order accurate having two nodes.

But if we want higher order accuracy at a point, let say 4^th order accuracy, then we will require more number of nodes to get this accuracy.

At each node, a taylor series is performed.

Taylor table makes it possible to find out any order accuracy of a derivative by taking required stencil points according to the formula.

For Asymmetric stencil ,

Order of approximation = (number of nodes in a stencil) – (order of derivative)

For Symmetric stencil

Order of approximation = (number of nodes in stencil) – (order of derivative) +1

Now, to find 4^th order accurate for second derivative by central differencing scheme.

If the node at which derivative function is to evaluate is considered as a central node and if equal number of nodes are taken from both the sides of a central node, then this method is said to be a central differencing and it is symmetric in nature.

For fourth order central differencing of second derivative

We have formula

Order of approximation = (number of nodes in a stencil ) - (order of derivative) +1

4 = (number of nodes in a stencil ) - (2) +1

Number of node in a stencil = 5

Hence to achieve 4^th order accuracy at a point for second derivative by central differencing, ( A stencil of five nodes is required)

	f(x)	f'(x)	f''(x)	f'''(x)	f''''(x)	f'''''(x)
a*f(x - 2h)	a	-a(2h)	a(2h)^2/2	-a(2h)^3/6	a (2h)^4/24	- a(2h)^5/120
b*f(x - h)	b	-bh	b(h)^2/2	-b(h)^3/6	b(h)^4/24	- b(h)^5/120
c*f(x)	c	0	0	0	0	0
d*f(x +h)	d	dh	d(h)^2/2	d(h)^3/6	d(h)^4/24	d(h)^5/120
e*f(x + 2h)	e	e(2h)	e(2h)^2/2	e(2h)^3/6	e(2h)^4/24	e(2h)^5/120
	0	0	1	0	0	?

we observe that, h term exists with every coefficient, so taking h common in their respective column, the taylor table becomes

	f(x)	(h)*f'(x)	(h)^2*f''(x)	(h)^3*f'''(x)	(h)^4*f''''(x)	(h)^5*f'''''(x)
a*f(x - 2h)	a	-2a	(4/2)a	-(8/6)a	(16/24)a	-(32/120)a
b*f(x - h)	b	-b	b/2	-b/6	b/24	-b/120
c*f(x)	c	0	0	0	0	0
d*f(x +h)	d	d	d/2	d/6	d/24	d/120
e*f(x + 2h)	e	2e	4e/2	8e/6	(16/24)e	(32/120)e
	0	0	1	0	0	?

by adding the columns, we get a hypothesis equation

a*f(x-2h) + b*f(x-h) + c*f(x) + d*f(x+h) + e*f(x +2h) = h^2*f''(x)

dividing by (h)^2 on both sides, we get

f''(x) = (a*f(x-2h) + b*f(x-h) + c*f(x) + d*f(x+h) + e*f(x +2h))/(h^2)

solving the taylor table by matrix inversion method in matlab to get the coefficients

% The below equations are for Central_difference_fourth_order_accurate for
% five nodes stencil


% a + b + c + d + e = 0

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

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

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

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


% Solving by AX = B , X = inv(A)*B


A = [ 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];

B = [0;0;1;0;0];


X = inv(A)*B

ln command window

>> coefficients_Central

X =

   -0.0833
    1.3333
   -2.5000
    1.3333
   -0.0833

hence we get,

a = -0.0833

b = 1.3333

c = -2.5000

d = 1.3333

e = -0.0833

using the function in matlab to get fourth order approximation for second derivative by central differencing.

the given function is f(x) = exp(x)*cos(x)

function outC = Central_fourth_order_approximation_for_second_derivative(x,h)

% analytical derivative 
% f"(x) = -2*exp(x)*sin(x)
  
x = pi/3;  
analytical_derivative = -2*exp(x)*sin(x);          

% numerical derivative, fourth_order_approximation for second derivative

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

  
  Central_differencing = (a*exp(x-2*h)*cos(x-2*h) + b*exp(x-h)*cos(x-h) + c*exp(x)*cos(x) + d*exp(x+h)*cos(x+h) + e*exp(x+2*h)*cos(x+2*h))/(h^2);
 
  outC  =  abs(Central_differencing - analytical_derivative);

end
    
    
    
    

Here the name of a function should be the name of a file.

Now to find 4th order accuracy of second order derivative of a given function for skewed right side differencing

If the nodes are considered on the axis and the neighbouring node points are considered only from the right side of the axis then it is called as skewed right side differencing.

In this, the nodes are asymmetric in nature

For Asymmetric stencil

order of approximation = (number of stencil points) - (order of derivative)

4 = (number of stencil points) - 2

Number of stencil points = 6.

	f(x)	(h)*f'(x)	(h)^2*f''(x)	(h)^3*f'''(x)	(h)^4*f''''(x)	(h)^5*f'''''(x)	(h)^6*f''''''(x)
a*f(x)	a	0	0	0	0	0	0
b*f(x+h)	b	b	b/2	b/6	b/24	b/120	b/720
c*f(x+2h)	c	2c	(4/2)c	(8/6)c	(16/24)c	(32/120)c	(64/720)c
d*f(x +3h)	d	3d	(9/2)d	(27/6)d	(81/24)d	(243/120)d	(729/720)d
e*f(x +4h)	e	4e	(16/2)e	(64/6)e	(256/24)e	(1024/120)e	(4096/720)e
f*f(x+5h)	f	5f	(25/2)f	(125/6)f	(625/24)f	(3125/120)f	(15625/720)f
	0	0	1	0	0	0	?

adding all the columns to get a hypothesis equation

a*f(x) + b*f(x +h) + c *f(x +2h) + d*f(x +3h) + e*f(x +4h) +f*(x +5h) = h^2*f''(x)

dividing by h^2, we get,

f''(x) = (a*f(x) + b*f(x +h) + c *f(x +2h) + d*f(x +3h) + e*f(x +4h) +f*(x +5h))/(h^2)

solving the taylor table of skewed right differencing by matrix inversion method in matlab to get the coefficients

% the below equations are of right skewed diiferencing
% a + b + c + d + e + f = 0
% 0 + b +2c + 3d + 4e +5f = 0
% 0 + b/2 + (4/2)c + (9/2)d + (16/2)e + (25/2)f = 1
% 0 + b/6 +(8/6)c + (27/6)d +(64/6)e + (125/6)f = 0
% 0 + b/24 + (16/24)c + (81/24)d + (256/24)e + (625/24)f = 0
% 0 + b/120 +(32/120)c +(243/120)d + (1024/120)e + (3125/120)f = 0


A = [1 1 1 1 1 1; 0 1 2 3 4 5; 0 1/2 4/2 9/2 16/2 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];

T = inv(A)*(B)

in command window

>> coefficients_Right_side

T =

    3.7500
  -12.8333
   17.8333
  -13.0000
    5.0833
   -0.8333

a = 3.7500

b = -12.8333

c = 17.8333

d = -13.0000

e = 5.0833

f = -0.8333

matlab code using function

function outR = rightskewed_Fourth_order_approximation_for_second_derivative(x,h)

% analytical derivative 
% f"(x) = -2*e^x*sin(x)
  
x = pi/3;  
analytical_derivative = -2*exp(x)*sin(x);

% numerical derivative, Skewed_fourth_order_approximation
  a = 3.7500;
  b  = -12.8333;
  c = 17.8333;
  d  = -13.0000;
  e =  5.0833;
  f = -0.8333;
   
% Hypothesis function
  skewed_right_side_forward_differencing = (a*exp(x)*cos(x) + b*exp(x+h)*cos(x+h) + c*exp(x+2*h)*cos(x+2*h) +d*exp(x+3*h)*cos(x+3*h)+ e*exp(x+4*h)*cos(x+4*h)+f*exp(x+5*h)*cos(x+5*h))/(h^2);
 
  outR =  abs(skewed_right_side_forward_differencing - analytical_derivative);

end

Now, to find 4th order accuracy of second order derivative of a given function for skewed left side differencing

If the nodes are considered on the axis and the neighbouring node points are considered only from the left side of the axis then it is called as skewed left side differencing.

In this, the nodes are asymmetric in nature

For Asymmetric stencil

order of approximation = (number of stencil points) - (order of derivative)

4 = (number of stencil points) - 2

Number of stencil points = 6.

	f(x)	h*f'(x)	(h)^2*f''(x)	(h)^3*f'''(x)	h^4*f''''(x)	(h)^5*f'''''(x)	(h)^6*f''''''(x)
a*f(x)	a	0	0	0	0	0	0
b*f(x-h)	b	-b	b/2	-b/6	b/24	- b/120	b/720
c*f(x-2h)	c	-2c	(4/2)c	-(8/6)c	(16/24)c	-(32/120)c	(64/720)c
d*f(x-3h)	d	-3d	(9/2)d	-(27/6)d	(81/24)d	-(243/120)d	(729/720)d
e*f(x-4h)	e	-4e	(16/2)e	-(64/6)e	(256/24)e	-(1024/120)e	(4096/720)e
f*f(x-5h)	f	-5f	(25/2)f	-(125/6)f	(625/24)f	-(3125/120)f	(15625/720)f
	0	0	1	0	0	0	?

we get a hypothesis equation as

f''(x) = (a*f(x) + b*f(x -h) + c*(x - 2h) +d*(x -3h) + e*f(x -4h) + f*f(x -5h))/(h^2)

% Coefficent for skewed left differencing for fourth order approximation of
% second derivative

% a + b + c + d + e + f = 0
% 0 - b -2c -3d -4e -5f = 0
% 0 + b/2 + (4/2)c + (9/2)d + (16/2)e + (25/2)f = 1
% 0 - b/6 - (8/6)c - (27/6)d - (64/6)e - (125/6)f = 0
% 0 + b/24 + (16/24)c + (81/24)d + (256/24)e + (625/24)f = 0
% 0 - b/120 -(32/120)c -(243/120)d - (1024/120)e - (3125/120)f = 0

 %Solving by A*X = B , X = inv(A)*B
 
 
A = [1 1 1 1 1 1; 0 -1 -2 -3 -4 -5; 0 1/2 2 9/2 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

In command

>> coefficients_Left_side

X =

    3.7500
  -12.8333
   17.8333
  -13.0000
    5.0833
   -0.8333

>> 

hence, we get

a = 3.7500

b = -12.8333

c = 17.8333

d = -13.0000

e = 5.0833

f = -0.8333

matlab code using function for left side differencing

function outL = leftskewed_fourth_order_approximation_for_second_derivative(x,h)

% analytical derivative 
% f"(x) = -2*e^x*sin(x)
  
x = pi/3;  
analytical_derivative = -2*exp(x)*sin(x);

% numerical derivative, fourth_order_approximation
 a = 3.7500;
 b = -12.8333;
 c = 17.8333;
 d = -13.0000;
 e =  5.0833;
 f = -0.8333;
 

 skewed_left_side_differencing = (a*exp(x)*cos(x) + b*exp(x-h)*cos(x-h) + c*exp(x-2*h)*cos(x-2*h) + d*exp(x-3*h)*cos(x-3*h)+ e*exp(x-4*h)*cos(x-4*h)+ f*exp(x-5*h)*cos(x-5*h))/(h^2);
 
  outL =  abs(skewed_left_side_differencing - analytical_derivative);

end

Main program

clear 
close all
clc

    x = pi/3;
    h = linspace(pi/4000,pi/40,20);
   
    for i = 1:length(h)
    
    
     outC(i) =  Central_fourth_order_approximation_for_second_derivative(x,h(i));
    
     outR(i) =  rightskewed_Fourth_order_approximation_for_second_derivative(x,h(i));

     outL(i)  =  leftskewed_fourth_order_approximation_for_second_derivative(x,h(i));
    
     

    end
    
      loglog(h,outC,'b')
      hold on
      loglog(h,outR,'g')
      hold on
      loglog(h,outL,'r')
      hold on
     
      legend('central difference error','skewed right error','skewed left error','Location', 'northwest')
      xlabel('h')
      ylabel('error')

we observe in the (image1){image 1 attached} that the skewed right differencing error is higher as compared to error of central and skewed left differencing.

To make it more clear if we increase the range of h (i.e pi/4000 to pi/4) and counts only ten values in between the range, then we can clearly observe the difference in the errors in

(image 2){image attached}

Conclusion

- In central differencing, due to treating a node at centre, the mesh size becomes compact

- Also the error in central differencing is low as compared to the error of the other differencing schemes

- This results in fine mesh. So it can be concluded that the central differencing is better for finer mesh.

- It is also applicable in the problems which contain diffusion phenomena because diffusion occur in every direction

- skewed type schemes can be applied for propogation/Marching problems as the flow has one direction.