Evaluation of fourth order approximation of second derivative of a function using Taylor table method following central and skewed differencing scheme

Objective :

Deriving the 4th order approximations of the second order derivative using central differencing, skewed right sided difference and skewed left sided difference scheme. Also, writing a MATLAB program which uses the developed approximations to evaluate the second derivative of the given function and to make comparison between different schemes.

Description :

The derivatives of the function, which are difficult to evaluate analytically, are found using numerical schemes. It also provides a benefit to evaluate the derivatives computationally, since it convert continous form of equation into algebraic form. Though the first order derivative can be found by doing a little manipulation in taylor's series, but for higher order approximations, we need to use Taylors table method. This method is based on assuming a hypothetical equation based on the order of approximation and collecting the coefficients, obtained after expanding each term using taylor series, in a tabular form.

Derivation of 4th order approximation of second derivative :

(1) Central differencing -

Since the approximation is 4th order, we need to take information from 4 extra points. Also, it is central difference scheme, so we take two points on the right hand side and two points on the left hand side. The numeric stencil will be as shown below :

Assuming hypothetical equation for second derivative:

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

Expanding each term using Taylor series and collecting coefficients:

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

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

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

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

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

We have five equations and five variables. If we take further equations, it will lead to overspecifying the system which results in abstract results. Though to check the order of truncation error or the truncation error itself, we need to take few more equations.

Writing above equations in matrix form :

`[[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]][[a], [b], [c], [d], [e]]=[[0],[0],[1],[0],[0]]`

Above equation can be written as:

`AX=B`

`X=A^-1B`

Solving above equation, which can be done using MATLAB, we get the values of the coefficients:

a = -0.0833

b = 1.3333

c = -2.5

d = 1.3333

e = -0.0833

putting these value in assumed equation,

`Deltax^2(del^2y)/(delx^2)=-0.0833f(i-2)+1.3333f(i-1)-2.5f(i)+1.3333f(i+1)-0.0833f(i+2)`

`(del^2y)/(delx^2)=1/(Deltax^2)(-0.0833f(i-2)+1.3333f(i-1)-2.5f(i)+1.3333f(i+1)-0.0833f(i+2))`

Above is the 4th order central differencing approximation of second derivative

(2) Skewed right sided differencing -

As the numerical scheme followed is skewed towards right side of the original point; numeric stencil will be on the right side i.e. the extra points from where the information is taken will be on the right side as shown below :

Assuming hypothetical equation for second derivative:

`(Deltax)^2(del^2y)/(delx^2) = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)`

According to Taylor's table method, each term is expanded and the coefficents are tabulated below:

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

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

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

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

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

Converting above equations in matrix form :

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

Above equation can be written as:

`AX=B`

`X=A^-1B`

Solving above equation, which can be done using MATLAB, we get the values of the coefficients:

a = 2.9167

b = -8.6667

c = 9.5

d = -4.6667

e = 0.9167

putting these value in assumed equation,

`Deltax^2(del^2y)/(delx^2)=2.9167f(i)-8.6667f(i+1)+9.5f(i+2)-4.6667f(i+3)+0.9167f(i+4)`

`(del^2y)/(delx^2)=1/(Deltax^2)(2.9167f(i)-8.6667f(i+1)+9.5f(i+2)-4.6667f(i+3)+0.9167f(i+4))`

Above is the 4th order skewed right sided differencing approximation of second derivative.

(3) Skewed left sided differencing -

Numeric stencil:

Assuming hypothetical equation for second derivative:

`(Deltax)^2(del^2y)/(delx^2) = af(i-4)+bf(i-3)+cf(i-2)+df(i-1)+ef(i)`

According to Taylor's table method, each term is expanded and the coefficents are tabulated below:

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

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

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

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

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

Converting above equations in matrix form :

`[[1,1,1,1,1],[-4, -3,-2,-1,0],[8,9/2,2,1/2,0],[-32/3,-9/2,-4/3,-1/6,0],[32/3,27/8,2/3,1/24,0]][[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

Above equation can be written as:

`AX=B`

`X=A^-1B`

Solving above equation, which can be done using MATLAB, we get the values of the coefficients:

a = 0.9167

b = -4.6667

c = 9.5

d = -8.6667

e = 2.9167

putting these value in assumed equation,

`Deltax^2(del^2y)/(delx^2)=0.9167f(i-4)-4.6667f(i-3)+9.5f(i-2)-8.6667f(i-1)+2.9167f(i)`

`(del^2y)/(delx^2)=1/(Deltax^2)(0.9167f(i-4)-4.6667f(i-3)+9.5f(i-2)-8.6667f(i-1)+2.9167f(i))`

Above is the 4th order skewed left sided differencing approximation of second derivative.

MATLAB Code :

Following is a program in Matlab to evaluate the second order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations that are derived above:

Main Program :

clear all
close all
clc

% Analytical function = exp(x)*cosx
% Second derivative = -2*exp(x)*sinx

x = pi/3;
dx = linspace(pi/30,pi/3000,100);

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

% Second derivative using numerical method

% central_differencing =  (1/dx^2)*[-0.0833f(i-2)+1.3333f(i-1)-2.5f(i)+1.3333f(i+1)-0.0833f(i+2)]

central_differencing =  (1./dx.^2).*[-0.0833*f(x-2*dx)+1.3333*f(x-dx)-2.5*f(x)+1.3333*f(x+dx)-0.0833*f(x+2*dx)];

% skewed_right_derivative = (1/dx^2)*[2.9167f(i)-8.6667f(i+1)+9.5f(i+2)-4.6667f(i+3)+0.9167f(i+4)]

skewed_right_derivative = (1./dx.^2).*[2.9167*f(x)-8.6667*f(x+dx)+9.5*f(x+2*dx)-4.6667*f(x+3*dx)+0.9167*f(x+4*dx)];

% skewed_left_derivative = (1/dx^2)*[0.9167f(i-4)-4.6667f(i-3)+9.5f(i-2)-8.6667f(i-1)+2.9167f(i)]

skewed_left_derivative = (1./dx.^2).*[0.9167*f(x-4*dx)-4.6667*f(x-3*dx)+9.5*f(x-2*dx)-8.6667*f(x-dx)+2.9167*f(x)];

% Error calculation

error_central = abs(central_differencing-second_derivative);
error_right = abs(skewed_right_derivative-second_derivative);
error_left = abs(skewed_left_derivative-second_derivative);

% Comparison Result

fprintf('Second order derivative of the function f(x) = exp(x)*cos(x) at x=pi/3 and dx = pi/30 is:\n')
fprintf('Analytical derivative = %d\n', second_derivative)
fprintf('Central differencing = %d\n', central_differencing(1))
fprintf('Skewed right differencing = %d\n', skewed_right_derivative(1))
fprintf('Skewed left differencing = %d\n', skewed_left_derivative(1))

% Plotting result

figure(1)
hold on
plot(dx,error_central,'color','r')
plot(dx,error_right,'color','b')
plot(dx,error_left,'color','g')
xlabel('Grid Size dx')
ylabel('Error')
title({'Error in 4th Order Approximation of', 'Second Derivative'})
legend('Central differencing','Skewed rightside differencing','Skewed leftside differencing','Location','NorthWest')

Function :

function [y] = f(i)

y = exp(i).*cos(i);

end

Results :

The above curve shows that the error in the results obtained by central differencing scheme do not vary largely with the increment in grid size whereas the error in the results obtained by skewed right side differencing and skewed left side differencing shows abrupt changes with the increment in grid size.

Why a Skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

From the above results, it is clear that the results obtained by central differencing technique are less fluctuating with grid size and are more reliable. But there might be situations where we need to evaluate the derivatives at either of the boundaries of the system. In such cases, we can only take information from the nearer points which are either on the left or right side of the original point; points on both side of the original point is not available. Here, skewed scheme is used. Central differencing cannot be used to evaluate derivatives of any order at the boundary of the system and only skewed scheme is helpful in this scenario.

Conclusion:

The 4th order approximations of the second order derivative using central differencing, skewed right sided difference and skewed left sided difference scheme are derived. The MATLAB program results, which uses the developed approximations to evaluate the second derivative of the given function, make comparison between different schemes. The importance of central differencing scheme and skewed scheme is shown by defining their area of application and reliability of results with increament in grid size.