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

Deraviation of 4th order approximations of the second-order derivative using various schemes:-

The given function is

`exp(x)*cos(x)`

The second derivative is :-

`-2*exp(x)*sin(x)`

1. Central Difference Scheme:-

As given, no. of Nodes (N) required for 4rth order approximation (P = 4) of 2nd order deravative (Q = 2) is :-

P = N - Q + 1

so, N = P + Q -1

so, N = 4 + 2 -1 = 5

For the central difference with 5 nodes,

the 4rth order approximation of 2nd order differential using central difference scheme with unknown coefficients(a,b,c,d,e) can be expressed as below,

`{del^2f}/{delx^2} approx f''(x) + Deltax^4f^6(x)*`(constants)` = {af(x-2dx) + bf(x-dx) + cf(x) + df(x+dx) + ef(x+2dx)}/{Deltax^2}`

Using Taylor's Table:-

	`f(x)`	`(Deltax)f'(x)`	`(Deltax)^2f''(x)`	`(Deltax)^3f'''(x)`	`(Deltax)^4f^{iv}(x)`	`(Deltax)^5f^{v}(x)`	`(Deltax)^6f^{vi}(x)`
`af(x-2dx)`	`a`	`-2a`	`2a`	`{-4a}/{3}`	`{2a}/{3}`	`{-4a}/{15}`	`{4a}/{45}`
`bf(x-dx)`	`b`	`-b`	`{b}/{2}`	`{-b}/{6}`	`{b}/{24}`	`{-b}/{120}`	`{b}/{720}`
`cf(x)`	`c`	`0`	`0`	`0`	`0`	`0`	`0`
`df(x+dx)`	`d`	`d`	`{d}/{2}`	`{d}/{6}`	`{d}/{24}`	`{d}/{120}`	`{d}/{720}`
`ef(x+2dx)`	`e`	`2e`	`2e`	`{4e}/{3}`	`{2e}/{3}`	`{4e}/{15}`	`{4e}/{45}`
	0	0	1	0	0	?	?

Solving for unknown coefficients,

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

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

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

`{-4a}/{3} - b/6 + d/6 + {4e}/{3} = 0`

`{2a}/{3} + b/24 + d/24 + 2e/3 = 0`

Using linsolve to solve the above linear equations:-

Thus the value of a is -0.0833, b is 1.3333, c is -2.5, d is 1.3333, e is -0.833

2. Skewed Right-Sided Difference:-

The 4rth order approximation of 2nd order differential using Skewed Right-Sided difference scheme with unknown coefficients(a,b,c,d,e,g) can be expressed as below,

`{del^2f}/{delx^2} approx f''(x) + Deltax^4f^6(x)*`(constants)` = {af(dx) + bf(x+dx) + cf(x+2dx) + df(x+3dx) + ef(x+4dx) + gf(x+5dx)}/{Deltax^2}`

Proof of order of accuracy according to the given rule (The order of accuracy P = Number of nodes in the stencil N - order of derivative Q )

So, As per the rule:- Order of accuracy P = N - Q

N (Number of Nodes or Points) = 6 [We have taken no. of points to be 6]

Q (Order of derivative) = 2 [We are solving for 2nd order of derivative]

So, P (the order of accuracy) should be 6 - 2 = 4.

Also, the order of accuracy we are solving for is 4, so as per the statement, the order of accuracy has been proved.

Using Taylor's Table:-

	`f(x)`	`(Deltax)f'(x)`	`(Deltax)^2f''(x)`	`(Deltax)^3f'''(x)`	`(Deltax)^4f^(iv)(x)`	`(Deltax)^5f^(v)(x)`
`af(x)`	`a`	`0`	`0`	`0`	`0`	`0`
`bf(x+dx)`	`b`	`b`	`b/2`	`b/6`	`b/24`	`b/120`
`cf(x+2dx)`	`c`	`2c`	`4c/2`	`8c/6`	`16c/24`	`32c/120`
`df(x+3dx)`	`d`	`3d`	`9d/2`	`27d/6`	`81d/24`	`243d/120`
`ef(x+4dx)`	`e`	`4e`	`16e/2`	`64e/6`	`256e/24`	`1024e/120`
`gf(x+5dx)`	`g`	`5g`	`25g/2`	`125g/6`	`625g/24`	`3125g/120`
	0	0	1	0	0	0

Solving for unknown coefficients:-

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

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

`b/2 + 4c/2 + 9d/2 + 16e/2 + 25g/2 = 1`

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

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

`b/120 + 32c/120 + 243d/120 + 1024e/120 + 3125g/120 = 0`

Solving above linear Equations we get :-

The coefficients thus calculated are :-

a = 3.75, b = -12.8333, c = 17.8333, d = -13, e= 5.0833, g = -0.8333

3. Skewed Left-Sided Difference:-

The 4rth order approximation of 2nd order differential using Skewed Left-Sided difference scheme with unknown coefficients(a,b,c,d,e,g) can be expressed as below,

`{del^2f}/{delx^2} approx f''(x) + Deltax^4f^6(x)*`(constants)` = {af(dx) + bf(x-dx) + cf(x-2dx) + df(x-3dx) + ef(x-4dx) + gf(x-5dx)}/{Deltax^2}`

Proof of order of accuracy according to the given rule (The order of accuracy P = Number of nodes in the stencil N - order of derivative Q )

So, As per the rule:- Order of accuracy P = N - Q

N (Number of Nodes or Points) = 6 [We have taken no. of points to be 6]

Q (Order of derivative) = 2 [We are solving for 2nd order of derivative]

So, P (the order of accuracy) should be 6 - 2 = 4.

Also, the order of accuracy we are solving for is 4, so as per the statement, the order of accuracy has been proved.

Using Taylor's Table:-

	`f(x)`	`(Deltax)f'(x)`	`(Deltax)^2f''(x)`	`(Deltax)^3f'''(x)`	`(Deltax)^4f^(iv)(x)`	`(Deltax)^5f^(v)(x)`
`af(x)`	`a`	`0`	`0`	`0`	`0`	`0`
`bf(x-dx)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`
`cf(x-2dx)`	`c`	`-2c`	`4c/2`	`-8c/6`	`16c/24`	`-32c/120`
`df(x-3dx)`	`d`	`-3d`	`9d/2`	`-27d/6`	`81d/24`	`-243d/120`
`ef(x-4dx)`	`e`	`-4e`	`16e/2`	`-64e/6`	`256e/24`	`-1024e/120`
`gf(x-5dx)`	`g`	`-5g`	`25g/2`	`-125g/6`	`625g/24`	`-3125g/120`
	0	0	1	0	0	0

Solving for the unknown coefficients we get:-

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

`- b - 2c - 3d - 4e - 5g = 0`

`b/2 + 4c/2 + 9d/2 + 16e/2 + 25g/2 = 1`

`-b/6 -8c/6 -27d/6 - 64e/6 -125g/6 = 0`

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

`-b/120 -32c/120 -243d/120 -1024e/120 -3125g/120 = 0`

The coefficients thus calculated are :-

a = 3.75, b = -12.8333, c = 17.8333, d = -13, e= 5.0833, g = -0.8333

The solution using matlab code is described below:-

A few Helper functions are defined:-

Function absError

function outputArg1 = absError(inputArg1,inputArg2)
%Calculates the absolute error
% 
outputArg1 = log(abs(inputArg1-inputArg2));
end

Function analyticDerivative

function [outputArg1] = analyticDerivitive(x)
% This function calcultates the second derivative of function 
% exp(x)*cos(x)
% Parameter x represents the point on which derivative is to be calculated.

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

end

Function centralDifference

function outputArg1 = centralDifference(x, dx, fun, X, analytic_result)
%This function calcultates the second deravative by central difference
%scheme. The parameters :-
% x = point at which derivative is to be calculated.
% dx = step size or grid size
% fun = given function
% X = coefficients for central scheme is calculated
% analytic_result = the value of derivative calculated by analytical method
value = (X(1)*fun(x-2*dx) + X(2)*fun(x-dx) + X(3)*fun(x) + X(4)*fun(x+dx) + X(5)*fun(x+2*dx))/(dx^2);
outputArg1 = absError(value, analytic_result);
end

Function leftDifference

function outputArg1 = leftDifference(x,dx, fun, X, analytic_result)
%This function calcultates the second deravative by left difference
%scheme. The parameters :-
% x = point at which derivative is to be calculated.
% dx = step size or grid size
% fun = given function
% X = coefficients for left sided scheme is calculated
% analytic_result = the value of derivative calculated by analytical method
value = (X(1)*fun(x)+X(2)*fun(x-dx)+ X(3)*fun(x-2*dx) + X(4)*fun(x-3*dx) + X(5)*fun(x-4*dx) + X(6)*fun(x-5*dx))/(dx^2);
outputArg1 = absError(value,analytic_result);
end

Function rightDifference

function outputArg1 = rightDifference(x,dx, fun, X, analytic_result)
%This function calcultates the second deravative by right difference
%scheme. The parameters :-
% x = point at which derivative is to be calculated.
% dx = step size or grid size
% fun = given function
% X = coefficients for right sided scheme is calculated
% analytic_result = the value of derivative calculated by analytical method
value = (X(1)*fun(x)+X(2)*fun(x+dx)+ X(3)*fun(x+2*dx) + X(4)*fun(x+3*dx) + X(5)*fun(x+4*dx) + X(6)*fun(x+5*dx))/(dx^2);
outputArg1 = absError(value, analytic_result);
end

Function givenFunction

function outputArg1 = givenFunction(x)
%Calculate the value of given function at point x.
%   Outputs the value of given function exp(x)*cos(x) at given point.
% Parameter x represents the point on which value is to be calculated.
outputArg1 = exp(x)*cos(x);
end

The driver Script

close all
clear all
clc

% Point on which derivative is to be calculated
x = pi/4;

% The varying Grid Size
dx = linspace(pi/10,pi/250,1000);

% The derivative calculated analyticaly
analytic_result = analyticDerivitive(x);

% Coefficients for central, left sided, right sided difference are calculated based on the matrix
% formed
% Central Difference Scheme
A_central = [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];
B_central = [0;0;1;0;0];
X_central = linsolve(A_central,B_central);

% Left Sided Difference Scheme
A_left = [1,1,1,1,1,1;0,-1,-2,-3,-4,-5;0,1/2,4/2,9/2,16/2,25/2;0,-1/6,-8/6,-27/6,-64/6,-125/6;0,1/24,16/24,81/24,256/24,625/24;0,-1/120,-32/120,-243/120,-1024/120,-3125/120];
B_left = [0;0;1;0;0;0];
X_left = linsolve(A_left,B_left);

% Right Sided Difference Scheme
A_right = [1,1,1,1,1,1;0,1,2,3,4,5;0,1/2,4/2,9/2,16/2,25/2;0,1/6,8/6,27/6,64/6,125/6;0,1/24,16/24,81/24,256/24,625/24;0,1/120,32/120,243/120,1024/120,3125/120];
B_right = [0;0;1;0;0;0];
X_right = linsolve(A_right,B_right);


% The errors matrix is initialized 
errors = zeros(length(dx),3);
for i = 1:length(dx)
    errors(i,1) = centralDifference(x,dx(i), @givenFunction, X_central, analytic_result);
    errors(i,2) = leftDifference(x,dx(i), @givenFunction, X_left, analytic_result);
    errors(i,3) = rightDifference(x,dx(i), @givenFunction, X_right, analytic_result);
    
end

figure(1)
plot(dx,errors(:,1), 'linewidth', 1.5);
title("Comparing different Schemes with Grid Size")
hold on
plot(dx,errors(:,2), 'linewidth', 1.5);
hold on
plot(dx, errors(:,3), 'linewidth', 1.5);
hold on
xlabel("Grid Size")
ylabel("Log(Error)")

legend('Central Difference Scheme', 'Left Difference Scheme', 'Right Difference Scheme')

The Error plots are below:-

Conclusions:-

The Central Difference method having less error even at higher grid sizes depicts its reliablity in accuracy as compared to other two schemes. As, it takes care of information from left as well as right side of the point of interest. Since there is lack of information at the boundary of domain from atleast one side of point of interest, so central difference scheme wont work.

So At boundary nodes either left or right sided skewed schemes are preferred.

As can be seen from the plot above skewed schemes have relatively same perfomance at smaller grid size.