Deriving fourth order approximations using Taylor Table Method

OBJECTIVE: To derive the fourth order approximation of a second order derivative using central differencing scheme,skewed right difference scheme and skewed left difference scheme with the help of Taylor table and also to compare the obtained values with analytical derivative value for the given function.

DERIVING THE FOURTH ORDER APPROXIMATIONS:

1.Fourth order approximation of a second order derivative using central differencing scheme:

Numerical stencil for central differencing scheme

Let us consider the approximation

`(del^2u)/(delx^2)~~ (af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/(dx^2)+O(dx^4)` ---------equation (1)

`af(i-2) = af(i)-(af\'(i)*(2dx))/(1!)+(af\'\'(i)*(2dx)^2)/(2!)-(af\'\'\'(i).(2dx)^3)/(3!)+(af\'\'\'\'(i).(2dx)^4)/(4!)-----`

expand the above terms till fourth order

similarly expand all the terms present in eq (1). Now these terms are represented in the form of a table to solve for the coefficients of eq (1) and this can be solved using Taylor Table Method.

stencil/taylor series terms	`f(i)`	`dx*f\'(i)`	`dx^2*f\'\'(i)`	`dx^3*f\'\'\'(i)`	`dx^4*f\'\'\'\'(i)`
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
sum of the coefficients	0	0	1	0	0

From the above table when we sum the derivative terms we get the following set of equations

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

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

`2a+b/2+0+d/2+2e = 1`

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

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

the above set of linear equations have to be solved using matrix inversion method and hence are to be represented in the form of matrix i.e.

AX = B, then we get

`[[1,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

X = inv(A)*B then we have

from the above values we have a = -0.0833, b = 1.3333, c = -2.5000, d = 1.3333, e = -0.0833

substituting in eq (1) yields to our approximation

`((del^2u)/(delx^2))~~((-0.0833)(f(i-2))+(1.3333)(f(i-1))+(-2.5000)(f(i))+(1.3333)(f(i+1))+(-0.0833)(f(i+2)))/(dx^2)+O(dx^4)`

2.Fourth order approximation using skewed right difference scheme:

Numerical stencil for skewed diiference scheme

Let us consider the approximation

`(del^2u)/(delx^2)~~(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4))/(dx^2)+O(dx^4)`-------equation (2)

expanding each of the terms and putting them in the taylor table then we have our table as shown below:

stencil/taylor series terms	`f(i)`	`dx*f\'(i)`	`dx^2*f\'\'(i)`	`dx^3*f\'\'\'(i)`	`dx^4*f\'\'\'\'(i)`
af(i)	a	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24
cf(i+2)	c	2c	4c/2	8c/6	16c/24
df(i+3)	d	3d	9d/2	27d/6	81d/24
ef(i+4)	e	4e	16e/2	64e/6	256e/24
sum of coefficients	0	0	1	0	0

adding the derivative terms we get the equations as:

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

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

`0+b/2+2c+9/2d+8e=1`

`0+b/6+8/6c+27/6d+64/6e=0`

`0+b/24+16/24c+81/24d+256/24e=0`

the above set of linear equations have to be solved using matrix inversion method and hence are to be represented in the form of matrix i.e.

AX = B, then we get

`[[1,1,1,1,1],[0,1,2,3,4],[0,1/2,2,9/2,8],[0,1/6,8/6,27/6,64/6],[0,1/24,1/24,81/24,256/24]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

X = inv(A)*B then we have

from the above we have a=2.9167, b = -8.6667, c = 9.5000, d = -4.6667, e = 0.9167. Substituting these values in eq (2) we have our approximation

`(del^2u)/(delx^2)~~(2.9167f(i)-8.6667f(i+1)+9.5000f(i+2)-4.6667f(i+3)+0.9167f(i+4))/(dx^2)+O(dx^4)`

3.Fourth order approximation using skewed left difference scheme:

Numerical stencil for skewed left differencing scheme

Let us consider the approximation

`(del^2u)/(delx^2)~~(af(i-4)+bf(i-3)+cf(i-2)+df(i-1)+ef(i))/(dx^2)+O(dx^4)` -------equation (3)

expanding each of the terms and putting them in the taylor table then we have our table as shown below:

stencil/taylor series terms	`f(i)`	`dx*f\'(i)`	`dx^2*f\'\'(i)`	`dx^3*f\'\'\'(i)`	`dx^4*f\'\'\'\'(i)`
af(i-4)	a	-4a	8a	-64a/6	-256a/24
bf(i-3)	b	-3b	9b/2	-27b/6	81b/24
cf(i-2)	c	-2c	2c	-8c/6	16c/24
df(i-1)	d	-d	d/2	-d/6	d/24
ef(i)	e	0	0	0	0
sum of coefficients	0	0	1	0	0

adding the derivative terms we get the equations as:

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

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

`8a+9/2b+2c+d/2+0=0`

`-64/6a-27/6b-8/6c-d/6+0=0`

`256/24a+81/24b+16/24c+d/24+0=0`

the above set of linear equations have to be solved using matrix inversion method and hence are to be represented in the form of matrix i.e.

AX = B, then we get

`[[1,1,1,1,1],[-4,-3,-2,-1,0],[8,9/2,2,1/2,0],[-64/6,-27/6,-8/6,-1/6,0],[256/24,81/24,16/24,1/24,0]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

X = inv(A)*B then we have

from the above we have a=0.9167, b = -4.6667, c = 9.5000, d = -8.6667, e = 2.9167. Substituting these values in eq (3) we have our approximation

`(del^2u)/(delx^2)~~(0.9167f(i-4)-4.6667f(i-3)+9.5000f(i-2)-4.6667f(i-1)+0.9167f(i))/(dx^2)+O(dx^4)`

4. Approximation using Forward differencing scheme:

The approximation for forward differencing scheme is

`(d^2f(x))/(dx^2)=2(f(x+a)-f(x)-f\'(x)*dx)/(dx^2)`

5. Approximation using Backward differencing scheme:

The approximation for backward differencing scheme is

`(d^2f(x))/(dx^2)=2(f(x-a)-f(x)+f\'(x)*dx)/(dx^2)`

PROGRAMS:

I. Function code for central differencing scheme:

% Defining the function for central differencing scheme

function error_central_diff = central_difference(x,dx)

 % f(x) = e^(x)*cosx is the given function
  
 % second order derivative for the given function
 
   second_analytical_derivative = -2*exp(x)*sin(x);
   
    % calculating the central differencing value
    
    central_differencing = (((-0.0833)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))+((1.3333)*((exp(x-dx))*(cos(x-dx))))-((2.5)*((exp(x))*(cos(x))))+((1.3333)*((exp(x+dx))*(cos(x+dx))))+((-0.0833)*((exp(x+(2*dx)))*(cos(x+(2*dx))))))/(dx^2);
    
    % calculating the absolute error
    
    error_central_diff = abs(central_differencing - second_analytical_derivative)
end

II.Function code for skewed right differencing scheme:

% Defining the function for skewed right differencing scheme

function error_skewed_right = skewed_right_difference(x,dx)

 % f(x) = e^(x)*cosx is the given function
  
 % second order derivative for the given function
 
 second_analytical_derivative = -2*exp(x)*sin(x);
 
 % Calculating the skewed right differencing value
   
    skewed_right_diff = (((2.9167)*((exp(x))*(cos(x))))-((8.6667)*((exp(x+dx))*(cos(x+dx))))+((9.5000)*((exp(x+(2*dx))*(cos(x+(2*dx))))))-((4.6667)*((exp(x+(3*dx)))*(cos(x+(3*dx)))))+((0.9167)*((exp(x+(4*dx)))*(cos(x+(4*dx))))))/(dx^2);
    
    % Calculating the absolute error
    
    error_skewed_right = abs(skewed_right_diff-second_analytical_derivative)
end

III. Function code for skewed left differencing scheme:

% Defining the function for skewed left differencing scheme

function error_skewed_left = skewed_left_difference(x,dx)

% f(x) = e^(x)*cosx is the given function
  
 % second order derivative for the given function

   second_analytical_derivative = -2*exp(x)*sin(x);
   
    % Calculating the skewed right differencing value
    
    skewed_left_diff = (((0.9167)*((exp(x-(4*dx))*(cos(x-(4*dx))))))-((4.6667)*((exp(x-(3*dx))*(cos(x-(3*dx))))))+((9.5000)*((exp(x-(2*dx))*(cos(x-(2*dx))))))-((8.6667)*((exp(x-dx))*(cos(x-dx))))+((2.9167)*((exp(x))*(cos(x)))))/(dx^2);
    
    % Calculating the absolute error
    
    error_skewed_left = abs(skewed_left_diff-second_analytical_derivative)
end

IV.Function code for forward differencing scheme:

% Defining the function for forward differencing scheme

function error_forward_diff = forward_difference(x,dx)

% f(x) = e^(x)*cosx is the given function

% first_analytical_derivative = exp(x)*(cos(x)-sin(x))
  
 % second order derivative for the given function
 
   second_analytical_derivative = -2*exp(x)*sin(x);
   
    % calculating the forward differencing value
   
   forward_differencing = (2*((exp(x+dx)*cos(x+dx))-(exp(x)*cos(x))-((exp(x)*(cos(x)-sin(x)))*dx)))/(dx^2);
   
    % calculating the absolute error
   
   error_forward_diff = abs(forward_differencing - second_analytical_derivative)
end

V. Function code for backward differencing scheme:

% Defining the function for backward differencing scheme

function error_backward_diff = backward_difference(x,dx)

% f(x) = e^(x)*cosx is the given function

% first_analytical_derivative = exp(x)*(cos(x)-sin(x))
  
 % second order derivative for the given function
 
   second_analytical_derivative = -2*exp(x)*sin(x);
   
   % calculating the backward differencing value
   
   backward_differencing = (2*((exp(x-dx)*cos(x-dx))-(exp(x)*cos(x))+((exp(x)*(cos(x)-sin(x)))*dx)))/(dx^2);
   
   % calculating the absolute error
   
   error_backward_diff = abs(backward_differencing - second_analytical_derivative)
   
end

VI. MAIN PROGRAM CODE:

clear all
close all
clc

% Declaring the value for \'x\'

x = pi/3;

% Declaring the variable \'dx\'

dx = linspace(pi/4,pi/40,10);

% To begin a for loop to calculate the values of different schemes

for i = 1:length(dx)   % declaring the variable \'i\'
    
    central_diff_value(i) = central_difference(x,dx(i));  % Calling the function central difference
    
    skewed_right_value(i) = skewed_right_difference(x,dx(i));  % Calling the function skewed right difference
    
    skewed_left_value(i) = skewed_left_difference(x,dx(i));     % Calling the function skewed left difference
    
    forward_diff_value(i) = forward_difference(x,dx(i));        % Calling the function forward difference
    
    backward_diff_value(i) = backward_difference(x,dx(i));       % Calling the function backward difference
    
end


% Plotting the graph using log-log plot

figure(1)

loglog(dx,central_diff_value,\'b\',\'linewidth\',3)
hold on
loglog(dx,skewed_right_value,\'r\',\'linewidth\',3)
hold on
loglog(dx,skewed_left_value,\'g\',\'linewidth\',3)
hold on
loglog(dx,forward_diff_value,\'c\',\'linewidth\',3)
hold on
loglog(dx,backward_diff_value,\'k\',\'linewidth\',3)

xlabel(\'dx\',\'FontSize\',20)
ylabel(\'error\',\'FontSize\',20)

legend(\'CENTRAL DIFF ERROR\',\'SKEWED RIGHT ERROR\',\'SKEWED LEFT ERROR\',\'FORWARD DIFF ERROR\',\'BACKWARD DIFF ERROR\')


    

OUTPUT:

ERRORS OUTPUT WHEN dx = pi/4

ERRORS OUTPUT WHEN dx = pi/40

From the above outputs we can observe that for small value of \'dx\', CDS gives the least value of error and FDS gives the highest value of error.

REPRESENTATION OF LOG-LOG PLOT:

CONCLUSION:

1. From the above plots we can say that Central differencing scheme gives the least error value when comapred to other differencing schemes. Central differencing scheme utilises the information available at next neighbouring points on both the sides from the point of consideration to give the approximation value. CDS results in giving higher order approximations and we know that as order of approximation increases the error decreases. The advantage of CDS is that its convergence rate is faster than all other differencing schemes.

2. Next scheme after CDS is Skewed Left Differencing Scheme and then comes Skewed Right Differencing Scheme in the context of least errors. Though initially Skewed Right Differencing Scheme is resulting in highest errors as per the plot but as \'dx\' value decreases it stands next to Skewed Left Differencing Scheme. Skewed differencing schemes utilises the information available at next neighbouring points on one side from the point of consideration either left or right to give the approximation value. This scheme is used in computing higher order approximations. When compared with CDS, the error obtained with skewed schemes are relatively high.

3. Backward differencing and Forward differencing are next in the hierarchy of errors.These are used to compute the first order approximations. These schemes use the information available at next neighbouring point on one side from the point of consideration to give the approximation value and hence the error obtained with these schemes are relatively high.

Advantage of Skewed schemes over Central differencing scheme:

Irrespective of any order of approximation, CDS gives least error value and this has been proved in the above discussions. But when it comes to typical boundary problems as shown below we cannot depend on Central differencing scheme rather we have to go with Skewed Differencing Scheme.

Here in this problem if we have to compute X₀ value using Central differencing scheme then we have to know x_-1 value which is unknown. And hence when we use Skewed right differencing scheme we can compute the required value. In this aspect Skewed schemes are advantageous over the Central differencing scheme and thereby surpasses it.