Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

Derive the `4^(th)` order approximation of the second-order derivative using central difference scheme

let f be the function whose second derivative is approximated at `i^(th)` point using function values at (i-2), (i-1), i, (i+1), (i+2) points. let h be the distance between the consecutive points.

let `(d^2f)/dx^2=a*f_(i-2)+b*f_(i-1)+c*f_i+d*f_(i+1)+e*f_(i+2)+O(h^n)`

n is the order of the scheme and is to be determined.

By Taylor's series, `f_(i+1)=f_i+(h/(1!))f'_i+(h^2/(2!))f''_i+(h^3/(3!))f'''_i+.....`

Applying the above series to `f_(i-2), f_(i-1), and f_(i+2)` and summing them appropriately gives

`(d^2f)/dx^2=(a+b+c+d+e)*f_i+(-2ah-bh+dh+2eh)*f'_i+(2ah^2+(bh^2)/2+(dh^2)/2+2eh^2)*f''_i+......`

On equating the corresponding terms of derivatives appropriately it gives

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

`-2ah-bh+dh+2eh=0`

`2ah^2+(bh^2)/2+(dh^2)/2+2eh^2=1`

`(-8ah^3)/6-(bh^3)/6+(dh^3)/6+(8eh^3)/6=0`

`(16ah^4)/24+(bh^4)/24+(dh^4)/24+(16eh^4)/24=0`

Solving for the coefficients gives `a=-1/(12h^2), b=4/(3h^2), c=-5/(2h^2), d=4/(3h^2), and e=-1/(12h^2)`

substituting the above values gives n=4.

So `(d^2f)/(dx^2)=(-f_(i-2)/12+(4f_(i-1))/3-(5f_i)/2+(4f_(i+1))/3-f_(i+2)/12)/h^2+O(h^4)`

Thus the second derivative of a function using a central difference scheme (fourth-order) is derived using Taylor's table.

Following a similar approach to determine the second derivative using forward difference (right-sided) scheme.

let f be the function whose second derivative is approximated at `i^(th)` using function values at i, (i+1), (i+2), (i+3), (i+4), and (i+5) points. let h be the distance between the consecutive points.

let `(d^2f)/dx^2=a*f_i+b*f_(i+1)+c*f_(i+2)+d*f_(i+3)+e*f_(i+4)+g*f_(i+5)+O(h^n)`

a, b, c, d, e, and g are constants to be determined.

Expanding the function values at (i+1), (i+2), (i+3), (i+4), and (i+5) using Taylor's series and summing them.

`(d^2f)/dx^2=(a+b+c+d+e+g)*f_i+(bh+2ch+3dh+4eh+5gh)*f'_i+((bh^2)/2+(4ch^2)/2+(9dh^2)/2+(16eh^2)/2+(25gh^2)/2)*f''_i+......`

On equating the corresponding coefficients of derivatives on both sides, it gives

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

`bh+2ch+3dh+4eh+5gh=0`

`(bh^2)/2+(4ch^2)/2+(9dh^2)/2+(16eh^2)/2+(25gh^2)/2=1`

.

.

`(bh^5)/120+(32ch^5)/120+(243dh^5)/2+(1024eh^5)/120+(3125gh^5)/120=0`

The above 6 equations are solved to get the 6 constants.

They are `a=15/(4h^2), b=-77/(6h^2), c=107/(6h^2), d=-13/h^2, e=61/(12h^2), and g=-5/(6h^2)`

Substituting the above values gives n=4.

So `(d^2f)/dx^2=((15f_i)/4-(77f_(i+1))/6+(107f_(i+2))/6-13f_(i+3)+(61f_(i+4))/12-(5f_(i+5))/6)/h^2+O(h^4)`

Thus fourth-order right-sided difference scheme is derived.

Following a similar approach, a fourth-order left-sided difference scheme is derived.

So `(d^2f)/dx^2=((15f_i)/4-(77f_(i-1))/6+(107f_(i-2))/6-13f_(i-3)+(61f_(i-4))/12-(5f_(i-5))/6)/h^2+O(h^4)`

A function f for `e^x*cos(x)` is created to calculate its values at different input values. It is as follows.

function y=f(x)

    y=exp(x)*cos(x);

end

It is then used to calculate the second derivative for different dx values using various numerical schemes.

The MATLAB code is as follows.

close all
clear all
clc

% f=exp(x)*cos(x)
% Analytical derivative, f'=exp(x)*exp(x)*(-2*sin(x))

% let x=pi/4 and dx varies from pi/1000 to pi/10

x=pi/3;
dx=linspace(pi/4000,pi/10,20);

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

% 1 row and 20 columns, each column for different value of dx
numerical_forward=zeros(1,20);   
forward_error=zeros(1,20);

numerical_central=zeros(1,20);
central_error=zeros(1,20);

numerical_backward=zeros(1,20);
backward_error=zeros(1,20);

for i=1:length(dx)
    
    numerical_forward(i)=((15/4)*f(x)-(77/6)*f(x+dx(i))+(107/6)*f(x+2*dx(i))-13*f(x+3*dx(i))+(61/12)*f(x+4*dx(i))-(5/6)*f(x+5*dx(i)))/(dx(i)^2);
    forward_error(i)=abs(exact_value-numerical_forward(i));
    
    numerical_central(i)=(-(1/12)*f(x-2*dx(i))+(4/3)*f(x-dx(i))-2.5*f(x)+(4/3)*f(x+dx(i))-(1/12)*f(x+2*dx(i)))/(dx(i)^2);
    central_error(i)=abs(exact_value-numerical_central(i));
    
    numerical_backward(i)=((15/4)*f(x)-(77/6)*f(x-dx(i))+(107/6)*f(x-2*dx(i))-13*f(x-3*dx(i))+(61/12)*f(x-4*dx(i))-(5/6)*f(x-5*dx(i)))/(dx(i)^2);
    backward_error(i)=abs(exact_value-numerical_backward(i));
    
end
    
plot(dx,forward_error,dx,central_error,dx,backward_error);
legend('Forward','Central','Backward');
xlabel('Grid-length(dx)');
ylabel('Error');
title('Plot of error for different numerical schemes vs dx');

The generated plot is as follows.

In general, central difference schemes have better accuracy (or higher-order) than single-sided schemes, visible from the above plot where the central difference has less error than other schemes at all values of dx. But at boundaries, a central difference scheme can't be employed.

Single-sided schemes perform better in cases where the flow (or information) flows in a single particular direction.