Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Aim: 4th Order Approximation of a Second-Order Derivative using Taylors Table Method by Incorporating Central, Right-Skewed and Left-Skewed Numerical Schemes
1. Central Difference Scheme
Calculating Number of Nodes required for each scheme using the formula:
Number of nodes, n = p+q-1
where, n = number of nodes
p = order of derivaive = 2
q = order of approximation = 4
n = 2+4-1 = 5
Therefore, for each Numerical Scheme we require 5 nodes
*----------*----------*----------*----------*
(i-2) (i-1) (i) (i+1) (i+2)
*-----dx----*
f''(x),i = a*f(x),i-2 + b*f(x),i-1 + c*f(x),i + d*f(x),i+1 + e*f(x),i+2--------------------------1
Numerical Approximation f''(x),i = (a*f(x),i-2 + b*f(x),i-1 + c*f(x),i + d*f(x),i+1 + e*f(x),i+2)/dx2
from Taylor Series Approximation, we get
a(f(x),i-2) = af(x),i - af'(x),i*(-2dx)/1! + af''(x),i(-2dx)2/2! - af'''(x),i(-2dx)3/3! + af''''(x),i(-2dx)4/4! + af'''''(x),i(-2dx)5/5! + af''''''(x),i(-2dx)6/6!
b(f(x),i-1) = bf(x),i - bf'(x),i*(-dx)/1! + bf''(x),i(-dx)2/2! - bf'''(x),i(-dx)3/3! + bf''''(x),i(-dx)4/4! + bf'''''(x),i(-dx)5/5! + bf''''''(x),i(-dx)6/6!
c(f(x),i) = cf(x),i
d(f(x),i+1) = bf(x),i + bf'(x),i*(dx)/1! + bf''(x),i(dx)2/2! + bf'''(x),i(dx)3/3! + bf''''(x),i(dx)4/4! + bf'''''(x),i(dx)5/5! + bf''''''(x),i(dx)6/6!
e(f(x),i+2) = ef(x),i + ef'(x),i*(2dx)/1! + ef''(x),i(2dx)2/2! + ef'''(x),i(2dx)3/3! + ef''''(x),i(2dx)4/4! + ef'''''(x),i(2dx)5/5! + ef''''''(x),i(2dx)6/6!
Using Taylor series approximation of each term of eq 1, we generate the Taylors Table
| f(x),i | f'(x),i*dx | f''(x),i*dx2/2! | f'''(x),i*dx3/3! | f''''(x),i*dx4/4! | f'''''(x),i*dx5/5! | f''''''(x),i*dx6/6! | |
|---|---|---|---|---|---|---|---|
| a*f(x),i-2 | a | -2a | 4/2a | -8/6a | 16/24a | -32/125a | 64/720a |
| b*f(x),i-1 | b | -b | b/2 | -b/6 | b/24 | -b/120 | b/720 |
| c*f(x),i | c | 0 | 0 | 0 | 0 | 0 | 0 |
| d*f(x),i+1 | d | d | d/2 | d/6 | d/24 | d/120 | d/720 |
| e*f(x),i+2 | e | 2e | 4e/2 | 8e/6 | 16/24e | 32/125e | 64/720e |
| 0 | 0 | 1 | 0 | 0 | 0 | 0 | |
Writing the tabulated in systems of Algebric Equations
a+b+c+d+e = 0
-2a-b+0+d+2e = 0
4/2a+b/2+0+d/2+4e/2 = 1
-8/6a -b/6+0+d/6+8e/6 = 0
16/24a+b/24+0+d/24+16/24e = 0
This can be written as AX=B
Where X = A-1 * B
Doing Calculations on Matlab we get
a = -0.0833
b = 1.3333
c = -2.5000
d = 1.3333
e = -0.0833
The approximate value of second order derivative with 4th order Central Difference Approximation can be obtained by
f''(x),i = (a*f(x),i-2 + b*f(x),i-1 + c*f(x),i + d*f(x),i+1 + e*f(x),i+2)/dx2
2. Right-Skewed Scheme
Calculating Number of Nodes required for each scheme using the formula:
Number of nodes, n = p+q
where, n = number of nodes
p = order of derivaive = 2
q = order of approximation = 4
n = 2+4 = 6
*----------*----------*----------*----------*----------*
(i) (i+1) (i+2) (i+3) (i+4) (i+5)
f''(x),i = a*f(x),i + b*f(x),i+1 + c*f(x),i+2 + d*f(x),i+3 + e*f(x),i+4
Following same steps as done in Central Difference Scheme we get
| f(x),i | f'(x),i*dx | f''(x),i*dx2/2! | f'''(x),i*dx3/3! | f''''(x),i*dx4/4! | f'''''(x),i*dx5/5! | f''''''(x),i*dx6/6! | |
|---|---|---|---|---|---|---|---|
| a*f(x),i | a | 0 | 0 | 0 | 0 | 0 | |
| b*f(x),i+1 | b | b | b/2 | b/6 | b/24 | b/120 | |
| c*f(x),i+2 | c | 2c | 4c/2 | 8c/6 | 16c/24 | 32c/120 | |
| d*f(x),i+3 | d | 3d | 9d/2 | 27d/6 | 814/24d | 243/120e | |
| e*f(x),i+4 | e | 4e | 16e/2 | 64e/6 | 250e/24 | 1024e/120 | |
| g*f(x),i+5 | g | 5g | 25g/2 | 125g/6 | 625g/24 | 3125g/120 | |
| 0 | 0 | 1 | 0 | 0 | 0 |
Writing the system of Aljebric Equations
a+b+c+d+e+f+g = 0
0+b+2c+3d+4e+5g = 0
0+b/2+4c/2+9d/2+16e/2+25g/2 = 1
0+b/6+8c/6+27d/6+64e/6+125g/6 = 0
0+b/24+16c/24+814d/24+250e/24+625g/24 = 0
0+b/120+32c/120+243e/120+1024e/120+3125g/120 = 0
From Gauss Elimination
A X = B
X = inv(A)*B
a = 3.75
b = -12.833
c = 17.833
d = -13.00
e = 5.083
g = -0.833
The approximate value of second order derivative with 4th order Central Difference Approximation can be obtained by
f''(x),i = (a*f(x),i-2 + b*f(x),i-1 + c*f(x),i + d*f(x),i+1 + e*f(x),i+2)/dx2
3. Left-Skewed Scheme
Calculating Number of Nodes required for each scheme using the formula:
Number of nodes, n = p+q
where, n = number of nodes
p = order of derivaive = 2
q = order of approximation = 4
n = 2+4 = 6
*----------*----------*----------*----------*----------*
(i-5) (i-4) (i-3) (i-2) (i-1) (i)
f''(x),i = a*f(x),i-4 + b*f(x),i-3 + c*f(x),i-2 + d*f(x),i-1 + e*f(x),i
Following same steps as done in Central Difference Scheme we get:
| f(x),i | f'(x),i*dx | f''(x),i*dx2/2! | f'''(x),i*dx3/3! | f''''(x),i*dx4/4! | f'''''(x),i*dx5/5! | f''''''(x),i*dx6/6! | |
|---|---|---|---|---|---|---|---|
| a*f(x),i-5 | a | -5a | 25a/2 | -125a/6 | 625a/24 | -3125a/120 | |
| b*f(x),i-4 | b | -4b | 16b/2 | -64b/6 | 256b/24 | -1024b/120 | |
| c*f(x),i-3 | c | -3c | 9c/2 | -27c/6 | 81c/24 | -243c/120 | |
| d*f(x),i-2 | d | -2d | 4d/2 | -8d/6 | 16d/24 | -32d/120 | |
| e*f(x),i-1 | e | -e | e/2 | -e/6 | e/24 | -e/120 | |
| g*f(x),i | g | 0 | 0 | 0 | 0 | 0 | |
| 0 | 0 | 1 | 0 | 0 | 0 |
Writing the system of Aljebric Equations
a+b+c+d+e+f+g = 0
-5a-4b-3c-2d-e+0 = 0
25a/2+16b/2+9c/2+4d/2+e/2+0 = 1
-125a/6-64b/6-27c/6-8d/6-e/6+0 = 0
625a/24+256b/24+81c/21+16d/24+e/24+0 = 0
-3125a/120-1024b/120-243c/120-32d/120-e/120+0 = 0
From Gauss Elimination
A X = B
X = inv(A)*B
a = -0.833
b = 5.0833
c = -13.000
d = 17.833
e = -12.833
g = 3.750
The approximate value of second order derivative with 4th order Central Difference Approximation can be obtained by
f''(x),i = (a*f(x),i-2 + b*f(x),i-1 + c*f(x),i + d*f(x),i+1 + e*f(x),i+2)/dx2
Matlab Code
% Calculating coefficients for Numerical Schemes
% Cenral Difference Scheme
A1 = [1 1 1 1 1 ; -2 -1 0 1 2; 4/2 1/2 0 1/2 4/2; -8/6 -1/6 0 1/6 8/6; 16/24 1/24 0 1/24 16/24];
B1 = [0 0 1 0 0]';
X1 = A1\B1; % faster than inv(A)
a1 = X1(1,1);
b1 = X1(2,1);
c1 = X1(3,1);
d1 = X1(4,1);
e1 = X1(5,1);
fprintf('Coefficients for Central Difference Scheme:\n')
disp(X1)
% Right Skewed Scheme
A2 = [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];
B2 = [0 0 1 0 0 0]';
X2 = A2\B2;
a2 = X2(1,1);
b2 = X2(2,1);
c2 = X2(3,1);
d2 = X2(4,1);
e2 = X2(5,1);
g2 = X2(6,1);
fprintf('Coefficients for Right Skewed Scheme\n')
disp(X2)
% Left Skewed Scheme
A3 = [1 1 1 1 1 1; -5 -4 -3 -2 -1 0; 25/2 16/2 9/2 4/2 1/2 0; -125/6 -64/6 -27/6 -8/6 -1/6 0; 625/24 256/24 81/24 16/24 1/24 0; -3125/120 -1024/120 -243/120 -32/120 -1/120 0];
B3 = [0 0 1 0 0 0]';
X3 = A3\B2;
a3 = X3(1,1);
b3 = X3(2,1);
c3 = X3(3,1);
d3 = X3(4,1);
e3 = X3(5,1);
g3 = X3(6,1);
fprintf('Coefficients for Left Skewed Scheme\n')
disp(X3)
Writing Functions for Each Approximation
1. Central Difference Scheme
function error_central_difference = central_difference_approx(x,dx)
% f(x) = exp(x)*cos(x)
% f''(x) = -2*exp(x)*sin(x)
analytical_derivative = -2*exp(x)*sin(x); % second order derivative
numerical_derivative = (((-0.0833)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))+((1.3333)*((exp(x-dx))*(cos(x-dx))))-((2.5000)*((exp(x))*(cos(x))))+((1.3333)*((exp(x+dx))*(cos(x+dx))))-((0.0833)*((exp(x+(2*dx)))*(cos(x+(2*dx))))))/(dx^2);
error_central_difference = abs((analytical_derivative - numerical_derivative));
2. Right Skewed Approximation Writing Functions for Each Approximation
function error_right_skewed = right_skewed_approx(x,dx)
% f(x) = exp(x)*cos(x)
% f''(x) = -2*exp(x)*sin(x)
analytical_derivative = -2*exp(x)*sin(x); % second order derivative
numerical_derivative = (((3.7500)*((exp(x))*(cos(x))))-((12.8333)*((exp(x+dx))*(cos(x+dx))))+((17.8333)*((exp(x+(2*dx)))*(cos(x+(2*dx)))))-((13.0000)*((exp(x+(3*dx)))*(cos(x+(3*dx)))))+((5.0833)*((exp(x+(4*dx)))*(cos(x+(4*dx)))))-((0.8333)*((exp(x+(5*dx)))*(cos(x+(5*dx))))))/(dx^2);
error_right_skewed = abs((analytical_derivative - numerical_derivative));
3 Left Skewed Approximation
function error_left_skewed = left_skewed_approx(x,dx)
% f(x) = exp(x)*cos(x)
% f''(x) = -2*exp(x)*sin(x)
analytical_derivative = -2*exp(x)*sin(x); % second order derivative
numerical_derivative = (((-0.8333)*((exp(x-(5*dx)))*(cos(x-(5*dx)))))+((5.0833)*((exp(x-(4*dx)))*(cos(x-(4*dx)))))-((13.0000)*((exp(x-(3*dx)))*(cos(x-(3*dx)))))+((17.8333)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))-((12.8333)*((exp(x-dx))*(cos(x-dx))))+((3.7500)*((exp(x))*(cos(x)))))/(dx^2);
error_left_skewed = abs((analytical_derivative - numerical_derivative));
Main Code
clear all
close all
clc
% main code
% plotting errors in each approximation with varying values of dx
x = pi/3;
dx = linspace(pi/4,pi/4000,20);
analytical_derivative = -2*exp(x)*sin(x);
% initialising
error_central_difference = zeros(size(dx));
error_right_skewed = zeros(size(dx));
error_left_skewed = zeros(size(dx));
for i = 1:length(dx)
error_central_difference(i) = central_difference_approx(x,dx(i));
error_right_skewed(i) = right_skewed_approx(x,dx(i));
error_left_skewed(i) = left_skewed_approx(x,dx(i));
end
% plotting
figure(1)
loglog(dx, error_central_difference, 'LineWidth', 2, 'Color', 'r', 'DisplayName', 'Central Differencing');
hold on
loglog(dx, error_right_skewed, 'LineWidth', 2, 'Color', 'b', 'DisplayName', 'Right Skewed Scheme');
loglog(dx, error_left_skewed, 'LineWidth', 2, 'Color', 'g', 'DisplayName', 'Left Skewed Scheme');
title('Error variation with decreasing dx in different Numerical Schemes')
xlabel('Grid size dx')
ylabel('Error in Approximation')
legend(Location='northwest')
hold off
Displaying Output
1. Coefficients
'2. Error Plot
CONCLUSION:
As seen from the plot, the central differencing scheme has very less error compared to other schemes. Therefore, central differencing is better than other schemes because it uses information from both the sides from the point of consideration to give approximative and accurate values.
As the left and right skewed difference methods uses data only from one side of point of interest it gives less approximative and accurate values about the point.
WHY A SKEWED SCHEME IS USEFUL:
The system which doesn’t contain sufficient information on both side of the point of interest. In such cases we cannot use central differencing scheme so skewed scheme is useful in such cases.
WHAT CAN A SKEWED SCHEME DO THAT A CD SCHEME CANNOT DO:
CDS cannot not be employed at the boundary nodes due to unavailability of right and left nodes at the same time. In such cases the skewed differencing scheme is used.