"Taylors Table 4th Order Approximation 2nd Order Derivative"

OBJECTIVE

• To derive the following 4th order approximations scheme of the second-order derivative:
  
  1. Central Difference
  2. Skewed Right-Sided Difference
  3. Skewed Left-Sided difference
  
  
• To write a program in MATLAB to evaluate the second-order derivative of the analytical function, i.e. `f(x)=e^x×cos(x)`, and compare with the 3 numerical approximations derived earlier.
  
  
• Using a plot, compare the absolute error with respect to above-mentioned schemes.
  
  
• Compare above-mentioned schemes with Forward and Backward Differencing approximation schemes.
  
  
• The usefulness of the skewed scheme over the Central Differencing scheme.

Taylor’s Expansion:

`f(x+∆x)= f(x)+∂f/∂x ∆x+(∂^2 f)/(∂x^2 ) (∆x)^2/(2!)+(∂^3 f)/(∂x^3 ) (∆x)^3/(3!)+⋯++(∂^n f)/(∂x^n ) (∆x)^n/(n!)`

1^st Order, 2^nd Order, or Higher Order Derivatives:

`∂f/∂x , (∂^2 f)/(∂x^2 ) ,……… , (∂^n f)/(∂x^n )`

Discretization: PDE to Algebraic Equation

For deriving the higher-order approximations, it is not as simple to just apply some algebraic operation and get the results but rather it becomes too complex to convert derivate into algebraic equations.

To ease out, we use a handy tool called as ‘Taylor Table’ to derive. The process starts with assuming a polynomial based on terms i.e. grid points which depends on an order of approximation and derivative. The coefficients of the polynomial are unknowns which are inserted into the Taylor table and evaluation done either through simple manipulation or matrix calculation.

4^th Order Central Difference Scheme for the Second Order Derivative

Let us consider 5 grid points including x as central grid point, i.e.

`(x-2Δx),(x- Δx),(x),(x+Δx),and (x+2Δx)`

Based on selected grid points, a polynomial is assumed shown below:

`(∂^2 f)/(∂x^2 )=a.f(x-2∆x)+b.f(x-∆x)+c.f(x)+d.f(x+∆x)+e.f(x+2∆x)`

So, on coupling together terms to form equations

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

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

4a + b + d +4e = 1*2!

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

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

Above algebraic equations written into matrix form,

On solving matrix,

On substituting,

`f^'' (x).(∆x)^2+((64a+b+d+64e))/(6!) f^(vi) (x).(∆x)^6`

Since, pairs (a, e) and (b, d) are equal. The sixth derivative doesn’t get eliminated.

On further equating with polynomial and dividing both sides with `(∆x)^2`,

`(-f(x-2∆x)+16.f(x-∆x)-30.f(x)+16.f(x+∆x)-f(x+2∆x))/(12(∆x)^2 ) = f^('') (x)+O(∆x)^4`

`O(∆x)^4` represents the truncation error of fourth order approximation.

Hence, we have derived the fourth-order central difference approximation scheme of the second-order derivative mentioned below:

`(∂^2 f)/(∂x^2 )≈(-f(x-2∆x)+16.f(x-∆x)-30.f(x)+16.f(x+∆x)-f(x+2∆x))/(12(∆x)^2 )`

4^th Order Skewed Right-Sided Difference Scheme for the Second Order Derivative

Let us consider 5 grid points including x as leftmost grid point, i.e.

`(x),(x+Δx),(x+2Δx),(x+3Δx),and (x+4Δx)`

Based on selected grid points, a polynomial is assumed shown below.

`(∂^2 f)/(∂x^2 )=a.f(x)+b.f(x+∆x)+c.f(x+2∆x)+d.f(x+3∆x)+e.f(x+4∆x)`

So, on coupling together terms to form equations

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

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

`b+4c+9d+16e=1×2!`

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

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

Above algebraic equations written into matrix form,

On solving matrix,

On Substtuing,

`f^('') (x).(∆x)^2+((b+64c+729d+4096e))/(6!) f^(vi) (x).(∆x)^6`

Since, pairs (a, e) and (b, d) are equal. The sixth derivative doesn’t get eliminated.

On further equating with polynomial and dividing both sides with `(∆x)^2,`

`(35.f(x)-104.f(x-∆x)+114.f(x-2∆x)-56.f(x-3∆x)+11.f(x-4∆x))/(12(∆x)^2 ) =f^('') (x)+O(∆x)^4`

`O(∆x)^4` represents the truncation error of fourth order approximation.

Hence, we have derived the fourth-order skewed left-sided difference approximation scheme of the second-order derivative mentioned below:

`(∂^2 f)/(∂x^2 )≈(35.f(x)-104.f(x-∆x)+114.f(x-2∆x)-56.f(x-3∆x)+11.f(x-4∆x))/(12(∆x)^2 )`

Function and it’s 2^nd Order Derivative

The function has been shown below for which 2^nd order derivative is to be evaluated at point

`f(x)=e^x×cos(x)`

`(∂^2 f)/(∂x^2) =-2×e^x×sin⁡(x)`

The exact 2^nd order derivative is obtained at ‘pi/3’. While, approximate values of 2^nd order derivative is obtained by different schemes mentioned below and compared with exact solution:

• 4^th Order Central Differencing Scheme
• 4^th Order Skewed Right-Sided Differencing Scheme
• 4^th Order Skewed Left-Sided Differencing Scheme

The computed absolute error is plotted based on discretization for different 4^th order scheme. To be much clearer about the different schemes, the consistency plot is obtained which helps us to understand which scheme is better w.r.t slope of respective schemes.

MATLAB CODE STRUCTURE

MATLAB CODE

Main Program:

%f(x) = exp(x)*cos(x) defined through function named as 'f_4order'
clear all
close all
clc
%input variable, x defined here
x = pi/3;
% discretization of function by the range of dx
dx = linspace( pi/15 , pi/400 , 50 );
% actual second derivative of function, f"(x) = -2*exp(x)*sin(x)
dydx_actual = - 2 * exp(x) * sin(x) ;
for m = 1:length(dx)
    % evaluating function at all grid points i.e. '9' grid points
    for j = 1:9
        gp(j) = f_4order(x+((j-5)*dx(m)));
    end
    % 4th order approoximation schemes to compute 2nd order derivative
    dydx_4order(m,:) = scheme(gp,dx(m));
end
% absolute error calculated for each of the above applied schemes
error = abs( dydx_actual - dydx_4order );
% absolute error plot
f_plot(error,dx);

Defined Functions:

1. 'f(x)' function: "f_4order"

function out=f_4order(x)
out = exp(x)*cos(x);
end

2. '4th Order Schemes' function: "scheme"

function out = scheme(gp,dx) 
%central differencing @ 4th order aproximation
dydx_central_4 = (- gp(3) + 16*gp(4) - 30*gp(5) + 16*gp(6) - gp(7)) / (12*(dx^2)) ;
%skewed right sided differencing @ 4th order aproximation
dydx_right_4 = (35*gp(5) - 104*gp(6) + 114*gp(7) - 56*gp(8) + 11*gp(9)) / (12*(dx^2)) ;    
%skewed right sided differencing @ 4th order aproximation
dydx_left_4 = (35*gp(5) - 104*gp(4) + 114*gp(3) - 56*gp(2) + 11*gp(1)) / (12*(dx^2)) ;   
out = [dydx_central_4 dydx_right_4 dydx_left_4];    
end

3. 'Plot' function: "f_plot"

function out = f_plot(error,dx) 
%line plot for 'absolute error' vs 'dx' to observe the variation for different schemes
figure (1)
plot( dx , error(:,1) , 'color' , 'b' , 'marker' , 'D' ) ;
hold on
plot( dx , error(:,2) , 'color' , 'k' , 'marker' , 's' ) ;
plot( dx , error(:,3) , 'color' , 'r' , 'marker' , 'o' ) ;
title( 'Variation of absolute error based on discretization @ different scheme' );
xlabel( 'dx' ) ;
ylabel( 'Absolute Error' ) ;
legend( 'Central Difference of 4th Order', 'Skewed Right Sided Difference of 4th Order', 'Skewed Left Sided Difference of 4th Order' )
%log plot for better understanding and observation of previous line plot 
figure (2)
loglog( dx , error(:,1) , 'color' , 'b' , 'marker' , 'D' ) ;
hold on
loglog( dx , error(:,2) , 'color' , 'k' , 'marker' , 's' ) ;
loglog( dx , error(:,3) , 'color' , 'r' , 'marker' , 'o' ) ;
title( 'Consistency of different scheme' );
xlabel('dx');
ylabel('Absolute Error');
legend('Central Difference of 4th Order', 'Skewed Right Sided Difference of 4th Order', 'Skewed Left Sided Difference of 4th Order' );
end

Comparison Plot between 4th Order Schemes for Second Order Derivative

The above plot shows the variation of absolute error for second order derivative based on discretization at different 4^th order schemes for the defined function. Through this plot, it is prominent to say that as grid space is decreased, the absolute error reduces for skewed right sided difference scheme of 4^th order. While, it’s not so obvious about the other two schemes irrespective of the absolute error being too low. The comment could be made on other two schemes by plotting the consistency plot. The below plot shows the consistency for different schemes based on discretization.

It can be observed that the central difference scheme would yield lower absolute error in comparison to other two schemes. It also approaches to zero absolute error at much higher rate. Thus, it must be preferred amongst the three 4^th order scheme.

Through the skewed left sided scheme plot, it can be seen that initially absolute error tends to reduce at much higher rate. As the ‘dx’ is reduced, the sudden increase in absolute error is observed and then reduction in absolute error at slower rate.

Comparison with Forward and Backward Difference Scheme

Forward Differencing Scheme

`∂f/∂x≈(f(x+∆x)-f(x))/(∆x)`

Backward Differencing Scheme

`∂f/∂x≈(f(x)-f(x-∆x))/(∆x)`

Skewed Scheme

The word skew means to be biased. Suppose, if the bulk of grid points are towards the right of the reference grid point, then it is said to be positively skewed or skewed right-sided scheme. But, if the bulk of grid points are towards the left of the reference grid point, then it is said to be negatively skewed or skewed left-sided scheme. Same thing could be observed and said regarding the forward and backward differencing scheme irrespective of only two grid points used to evaluate the first order derivative.

Skewed Right-Sided Scheme

`(∂^2 f)/(∂x^2 )≈(35.f(x)-104.f(x+∆x)+114.f(x+2∆x)-56.f(x+3∆x)+11.f(x+4∆x))/(12(∆x)^2 )`

Skewed Left-Sided Scheme

`(∂^2 f)/(∂x^2 )≈(35.f(x)-104.f(x-∆x)+114.f(x-2∆x)-56.f(x-3∆x)+11.f(x-4∆x))/(12(∆x)^2 )`

On observing above two schemes, we can say that the coefficient of function at particular pairs of grid point are same. These pairs can be (0, 0), (-1, +1), (-2, +2), (-3, +3), … and so on depending on the scheme.

Thus, any skewed scheme can be applied anywhere in the grid in the form as either a left-sided or right-sided scheme.

Benefit of Skewed Scheme over Central Differencing Scheme

In order to apply the approximation scheme at boundaries, only skewed schemes are applicable. While due to the unavailability of grid points on one of the sides of the boundary, the central difference scheme can’t be applied. Because CDS requires a number of grid points on both sides of the reference points depending on the order of the approximation scheme and order of derivative to be evaluated.