Taylor table method and Matlab coding

Objectives:

• To derive the 4th order approximations of second order derivatives using;

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

1) Fourth order approximation of second order derivative using central differencing :

Consider the numerical stencil to be around a point 'i' in the central differencing scheme.

Let the hypothesis function be,

`(delta^2f )/( deltax^2)` = a`*`f(i-2) + b`*`(f(i-1)+ c`*`f(i) +d`*`f(i+1) + e`*`f(i+2) -------------------(1)

From the Taylor's table, the values of the co efficients can be found.

Taylor's Table

	f(i)	Δ x*f’(i)	(Δx^2)*f’’(i)	(Δ x^3)*f’’’(i)	(Δ x^4)*`f^(iv)`(i)	(Δ x^5)*`f^(v)`(i)	(Δ x^6)*`f^(vi)`(i)
a*f(i-2)	a	-2*a	(4/2)*a	-(8/6)*a	(16/24)*a	-(32/120)*a	(64/720)*a
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	2*e	(4/2)*e	(8/6)*e	(16/24)*e	(32/120)*e	(64/720)*e
sum of the coefficients	0	0	1	0	0	?	?

The sum of coefficients of (Δx^2)*f''(i) alone is equal to one and all other sums are zero as only the second derivative of the function is to be found out. Also, we only need 5 equations as there are only 5 unknowns. so, the last 2 columns are ignored for now. Hence we have,

a + b + c + d + e = 0

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

2*a + b/2 + d/2 + 2*e = 1

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

(2/3)*a + b/24 + d/24 + (2/3)*e = o

Converting equation into matrix form,

A = [1 1 1 1 1; -2 -1 0 1 2; 2 0.5 0 0.5 2; -4/3 -1/6 0 1/6 4/3; 2/3 1/24 0 1/24 2/3]

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

X = AB

By solving the above equations we get,

a = -0.0833

b = 1.3333

c = -2.5000

d = 1.3333

e = -0.0833

Hence, a =e , b = d

Substituting these values in the entire taylor's table,

sum of column 2 = a + b + d + c + e = 0

similarly,

sum of column 3 = 0

sum of column 4 = 1

sum of column 5 = 0

sum of column 6 = 0

sum of column 7 = ((32/120)*a) -(b/120) + 0 +(b/120)-((32/120)*a) = 0

sum of column 8 =2*((64*a) +b)/720 = non zero number

Rewriting equation (1) after substituting for the a, b, c, d& e,

`(delta^2f )/( deltax^2)` = a*f(i-2) + b*(f(i-1)+ c*f(i) +d*f(i+1) + e*f(i+2) = ((Δx^2)*f’’(i)) + (2*((64*a)+b)/720)*`f^(vi)`(i)*(`Δ x^6`)

dividing LHS & RHS with Δx^2

(a*f(i-2) + b*(f(i-1)+ c*f(i) +d*f(i+1) + e*f(i+2))/Δx^2 = f’’(i)) + (2*((64*a) +b)/240)*`f^(vi)`(i)*(`Δ x^4`) ----------------------------------------------------------------(2)

We can rewrite the equation as,

(a*f(i-2) + b*f(i-1)+ c*f(i) +d*f(i+1) + e*f(i+2))/Δx^2 = f’’(i)) + O*(Δ x^4)

Therefore derived the fourth order approximation for second derivative in skewed right sided differencing scheme. Here, O*(Δ `x^4`) denotes the order of approximation and the truncation error.

2) Fourth order approximation of second order derivative using skewed right sided differencing :

Consider a numerical scheme as shown below,

Let the hypothesis function be,

`(delta^2f )/( deltax^2)` = a`*`f(i) + b`*`(f(i+1)+ c`*`f(i+2) +d`*`f(i+3) + e`*`f(i+4) + h`*`f(i+5)-------------------(1)

From the Taylor's table, the values of the co efficients can be found.

Taylor's Table

	f(i)	Δ x*f’(i)	(`Δx^2`)*f’’(i)	(`Δ x^3`)*f’’’(i)	(`Δ x^4`)*`f^(iv)`(i)	(`Δ x^5`)*`f^(v)`(i)	(`Δ x^6`)*`f^(vi)`(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	2*c	(4/2)*c	(8/6)*c	(16/24)*c	(32/120)*c	(64/720)*c
d*f(i+3)	d	3*d	(9/2)*d	(27/6)*d	(81/24)*d	(243/120)*d	(729/720)*d
e*f(i+4)	e	4*e	(16/2)*e	(64/6)*e	(256/24)*e	(1024/120)*e	(4096/720)*e
h*f(i+5)	h	5*h	(25/2)*h	(125/6)*h	(625/24)*h	(3125/120)*h	(15625/720)*h

As we have 6 unknowns now, so we need 6 equations to find the unknowns.

We have the sums of each columns as,

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

0 + b + 2*c + 3*d + 4*e + 5*h = 0

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

0 + b/6 + (4/3)*c + (27/6)*d + (64/6)*e +(125/6)*h = 0

0 + b/24 + (16/24)*c + (81/24)*d + (256/24)*e + (625/24)*h = 0

0 + b/120 + (`2^5`/120)*c + (`3^5`/120)*d + (`4^5`/120)*e + (`5^5`/120)*h = 0

Converting equation into matrix form,

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 4/3 27/6 64/6 125/6; 0 1/24 2/3 81/24 256/24 625/24; 0 1/120 (`2^5`/120) (`3^5/120`) (`4^5/120)` `(5^5/120)` ]

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

X = AB

By solving the above equations we get,

a = 3.75000

b = -12.83333

c = 17.83333

d = -13.0000

e = 5.08333

h = -0.83333

Substituting these values in the entire taylor's table,

sum of column 2 = a + b + d + c + e = 0

similarly,

sum of column 3 = 0

sum of column 4 = 1

sum of column 5 = 0

sum of column 6 = 0

sum of column 7 = 0

sum of column 8 =(b+ 64*c+ 729*d + 4096*e+ 15625*h )/720 = non zero number

Rewriting equation (1) after substituting for the a, b, c, d& e,

`(delta^2f )/( deltax^2)` =a`*`f(i) + b`*`(f(i+1)+ c`*`f(i+2) +d`*`f(i+3) + e`*`f(i+4) + h`*`f(i+5) = (`Δx^2`*f’’(i)) + ((b+ 64*c+ 729*d + 4096*e+ 15625*h )/720)*`f^(vi)`(i)*(`Δ x^6`)

dividing LHS & RHS with `Δx^2`

( a`*`f(i) + b`*`f(i+1)+ c`*`f(i+2) +d`*`f(i+3) + e`*`f(i+4) + h`*`f(i+5))/`Δx^2` = f’’(i)) + ((b+ 64*c+ 729*d + 4096*e+ 15625*h )/720)*`f^(vi)`(i)*(`Δ x^4`) ---------------(2)

Therefore derived the fourth order approximation for second derivative in skewed right sided differencing scheme as,

( a`*`f(i) + b`*`f(i+1)+ c`*`f(i+2) +d`*`f(i+3) + e`*`f(i+4) + h`*`f(i+5))/`Δx^2` = f’’(i)) + O*(`Δ x^4`)

3) Fourth order approximation of second order derivative using skewed left sided differencing :

Consider a numerical scheme as shown below,

Let the hypothesis function be,

`(delta^2f )/( deltax^2)` = a`*`f(i-5) + b`*`(f(i-4)+ c`*`f(i-3) +d`*`f(i-2) + e`*`f(i-1) + h`*`f(i)-------------------(1)

From the Taylor's table, the values of the co efficients can be found.

Taylor's Table

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

As we have 6 unknowns now, so we need 6 equations to find the unknowns.

We have the sums of each columns as,

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

0 - e - 2*d - 3*c - 4*b - 5*a = 0

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

0 - e/6 - (4/3)*d - (27/6)*c - (64/6)*b - (125/6)*a = 0

0 + e/24 + (16/24)*d + (81/24)*c + (256/24)*b + (625/24)*a = 0

0 - e/120 - (`2^5`/120)*d - (`3^5`/120)*c - (`4^5`/120)*b - (`5^5`/120)*a = 0

Converting equation into matrix form,

A = [1, 1, 1, 1, 1, 1; -5, -4, -3, -2, -1, 0; 25/2, 8, 9/2, 2, 1/2, 0; -125/6, -64/6, -27/6, -4/3, -1/6, 0; 625/24 256/24, 81/24, 2/3, 1/24, 0; -`(5^5/120)`, -(`4^5/120)`, -(`3^5/120`), -(`2^5`/120), -(1/120), 0]

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

X = AB

By solving the above equations we get,

a = -0.83333

b = 5.08333

c = -13.00000

d = 17.83333

e = -12.83333

h = 3.75000

Substituting these values in the entire taylor's table,

sum of column 2 = a + b + d + c + e = 0

similarly,

sum of column 3 = 0

sum of column 4 = 1

sum of column 5 = 0

sum of column 6 = 0

sum of column 7 = 0

sum of column 8 =(e + 64*d + 729*c + 4096*b + 15625*a )/720 = non zero number

Rewriting equation (1) after substituting for the a, b, c, d& e,

`(delta^2f )/( deltax^2)` =a`*`f(i-5) + b`*`(f(i-4)+ c`*`f(i-3) +d`*`f(i-2) + e`*`f(i-1) + h`*`f(i) = (`Δx^2`*f’’(i)) + ((e+ 64*d+ 729*c + 4096*b+ 15625*a )/720)*`f^(vi)`(i)*(`Δ x^6`)

dividing LHS & RHS with `Δx^2`

( a`*`f(i-5) + b`*`f(i-4)+ c`*`f(i-3) +d`*`f(i-2) + e`*`f(i-1) + h`*`f(i))/`Δx^2` = f’’(i)) + ((b+ 64*c+ 729*d + 4096*e+ 15625*h )/720)*`f^(vi)`(i)*(`Δ x^4`) --------------(2)

Therefore derived the fourth order approximation for second derivative in skewed left sided scheme as,

( a`*`f(i-5) + b`*`f(i-4)+ c`*`f(i-3) +d`*`f(i-2) + e`*`f(i-1) + h`*`f(i))`/``Δx^2` = f’’(i)) + O*(`Δ x^4`)

Thus derived the 4th order approximations of second order derivatives using different approximation techniques.

Program in Matlab to evaluate the second order derivative of the analytical function exp(x)*cos(x) and comparing it with the numerical approximations:

clear all
close all
clc

%values for x and dx
x = pi/4;
dx = linspace(pi/4,pi/4000,20);

for i = 1: length(dx)

 central_error(i) = central_diff_approximation(x, dx(i));  %Error by central differencing
 
 rightskewed_error(i) = rightskewed_approximation(x, dx(i));  %Error by skewed right sided differencing
 
 leftskewed_error(i) = leftskewed_approximation(x,dx(i));  %Error by skewed left sided differencing
 
 fds_error(i) = forward_approximation(x,dx(i));  %Error by forwad differencing
 
 bds_error(i) = backward_approximation(x,dx(i));  %Error by backward differencing
 
end  

%plotting the dx values against errors for each type of approximation
figure(1)
loglog(dx, central_error,'linewidth',2)
hold on
loglog(dx, rightskewed_error, 'r','linewidth',2)
hold on
loglog(dx, leftskewed_error,'g','linewidth',2)
hold on
loglog(dx, fds_error, 'y','linewidth',2)
hold on
loglog(dx, bds_error,'linewidth', 2)
grid
xlabel('dx')
ylabel('Absolute error')
legend('Central error', 'Right skewed error', 'Left skewed error', 'FDS error', 'BDS error')

Function for central difference scheme;

%Function for central differencing scheme 

function out = central_diff_approximation(x, dx)
  
  %f(x) = exp(x)*cos(x)
  %Analytical derivative = f'(x) = -2*exp(x)*sin(x) 
  
  analytical_derivative = -2*exp(x)*sin(x);
  
  % numerical derivative
  %fourth order central differencing approximation
  a = -0.0833;
  b = 1.3333;
  c = -2.5000;
  d = 1.3333;
  e = -0.0833;
  %central_differencing = (a*f(i-2) + b*f(i-1)+ c*f(i) +d*f(i+1) + e*f(i+2))/?x^2
  
 central_differencing = ((a*exp(x-(2*dx))*cos(x-(2*dx))) + (b*exp(x-dx)*cos(x-dx))+ (c*exp(x)*cos(x)) + (d*exp(x+dx)*cos(x+dx)) + (e*exp(x+(2*dx))*cos(x+(2*dx))))/(dx^2);
 
 %output
     out = abs(central_differencing - analytical_derivative)

endfunction

Function for right skewed difference scheme;

%Function for Skewed right sided differencing scheme 

function out = rightskewed_approximation(x, dx)
  
  %f(x) = exp(x)*cos(x)
  %Analytical derivative = f'(x) = -2*exp(x)*sin(x) 
  
  analytical_derivative = -2*exp(x)*sin(x);
  
  % numerical derivative
  %fourth order Skewed right sided differencing approximation
  a = 3.75000;
  b = -12.83333;
  c = 17.83333;
  d = -13.0000;
  e = 5.08333;
  h = -0.83333;
  %rightskewed_differencing = ( a*f(i) + b*f(i+1)+ c*f(i+2) +d*f(i+3) + e*f(i+4) + h*f(i+5))/?x^2
  
 rightskewed_differencing = ((a*exp(x)*cos(x)) + (b*exp(x+dx)*cos(x+dx))+ (c*exp(x+(2*dx))*cos(x+(2*dx))) + (d*exp(x+(3*dx))*cos(x+(3*dx))) + (e*exp(x+(4*dx))*cos(x+(4*dx))) + (h*exp(x+(5*dx))*cos(x+(5*dx))))/(dx^2);
 
 %output
     out = abs(rightskewed_differencing - analytical_derivative)

endfunction

Function for left skewed difference scheme;

%Function for Skewed left sided differencing scheme 

function out = leftskewed_approximation(x, dx)
  
  %f(x) = exp(x)*cos(x)
  %Analytical derivative = f'(x) = -2*exp(x)*sin(x) 
  
  analytical_derivative = -2*exp(x)*sin(x);
  
  % numerical derivative
  %fourth order Skewed left sided differencing approximation
  a = -0.83333;
  b = 5.08333;
  c = -13.00000;
  d = 17.83333;
  e = -12.83333;
  h = 3.75000;
  %leftskewed_differencing = ( h*f(i) + e*f(i-1)+ d*f(i-2) +c*f(i-3) + b*f(i-4) + a*f(i-5))/dx^2
  
 leftskewed_differencing = ((h*exp(x)*cos(x)) + (e*exp(x-dx)*cos(x-dx))+ (d*exp(x-(2*dx))*cos(x-(2*dx))) + (c*exp(x-(3*dx))*cos(x-(3*dx))) + (b*exp(x-(4*dx))*cos(x-(4*dx))) + (a*exp(x-(5*dx))*cos(x-(5*dx))))/(dx^2);
 
 %output
     out = abs(leftskewed_differencing - analytical_derivative)

endfunction

Function for forward difference scheme;

%Function for Forward differencing scheme 

function out = forward_approximation(x, dx)

   %f(x) = exp(x)*cos(x)
  %Analytical derivative = f'(x) = -2*exp(x)*sin(x) 
  
  analytical_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; 
   

    %output
     out = abs(forward_differencing - analytical_derivative)
  
  end

Function for backward difference scheme;

%Function for Backward differencing scheme 

function out = backward_approximation(x, dx)

   %f(x) = exp(x)*cos(x)
  %Analytical derivative = f'(x) = -2*exp(x)*sin(x) 
  
  analytical_derivative = -2*exp(x)*sin(x);
  
  % numerical derivative
  % backward_differencing = (f(x-dx) -  f(x))/dx
   
    backward_differencing = ((exp(x-dx)*cos(x-dx)) - (exp(x)*cos(x)))/dx; 
   

    %output
     out = abs(backward_differencing - analytical_derivative)
  
  end

Plot that compares the absolute error between the above mentioned schemes;

Conclusion

• From the above plot it is clear that the fourth order differencing approaches are very much more accurate than the first order difference approaches like FDS and BDS. For larger dx values, the errors caused by all the different difference approaches seems to be in the same range. But when dx is being reduced, fourth order approximations get closer to the exact value in an accelerated rate although the FDS and BDS shows no significant change in the errors. Hence proven that, higher order approximations are significantly more accurate than first order approximations.
• When comparing the fourth order approximations, we find that the skewed difference scheme gives the lowest errors than the central difference scheme. With the decrease in dx value, the errors of skewed difference schemes are also reducing. But for the central difference scheme, after some point, the error value is a constant even after further reduction in dx value.
• At times when there are no nodes on anyone of the either sides of the node we consider to discretize, we cannot use central differece scheme because it requires informations from nodes of either sides. In this type of situations, skewed difference schemes are useful because they require information only from any one of the sides of the node under consideration.