Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Skewed left-sided difference:
The order of accuracy P = Number of nodes in the stencil N -order of derivative Q
Number of nodes in the stencil N = The order of accuracy P+order of derivative Q
N = 4+2
N =6
Taylor Table for the skewed left side difference
| f(i) | `Deltaxf'(i)` | `Delta x^2f^2(i)` | `Delta x^3f^3(i)` | `Delta x^4f^4(i)` | `Delta x^5f^5(i)` | |
| af(i) | a | 0 | 0 | 0 | 0 | 0 |
| bf(i-1) | b | -b | `b/2` | `-b/6` | `b/24` | `-b/120` |
| cf(i-2) | c | -2c | `(4c)/2` | `(-8c)/6` | `(16c)/24` | `(-32c)/120` |
| df(i-3) | d | -3d | `(9d)/2` | `(-27d)/6` | `(81d)/24` | `(-243d)/120` |
| ef(i-4) | e | -4e | `(16e)/2` | `(-64e)/6` | `(256e)/24` | `(-1024e)/120` |
| gf(i-5) | g | -5g | `(25g)/2` | `(-125g)/6` | `(625g)/24` | `(-3125g)/120` |
| Sum | 0 | 0 | 1 | 0 | 0 | 0 |
Find the values of coefficient using MATLAB
M = [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];
N = [0;0;1;0;0;0];
Z = rats(inv(M)*N)
Z =
6×14 char array
' 15/4 '
' -77/6 '
' 107/6 '
' -13 '
' 61/12 '
' -5/6 'final equation = `(45f(i)-154f(i-1)+214f(i-2)-156f(i-3)+61f(i-4)-10f(i-5))/(12dx^2`
Skewed right-sided difference:
Number of nodes in the stencil N = The order of accuracy P+order of derivative Q
N = 4+2
N =6
| f(i) | `Deltaxf'(i)` | `Deltax^2f^2(i)` | `Delta x^3f^3(i)` | `Deltax^4f^4(i)` | `Deltax^5f^5(i)` | ||
| af(i) | a | 0 | 0 | 0 | 0 | 0 | |
| bf(i+1) | b | b | b/2 | b/6 | b/24 | b/120 | |
| cf(i+2) | c | 2c | 4c/2 | 8c/6 | 16c/24 | 32c/24 | |
| df(i+3) | d | 3d | 9d/2 | 27d/6 | 81d/24 | 243d/120 | |
| ef(i+4) | e | 4e | 16e/2 | 64e/2 | 256e/24 | 1024e/120 | |
| gf(i+5) | g | 5g | 25g/2 | 125g/6 | 625g/24 | 3125g/120 | |
| Sum | 0 | 0 | 1 | 0 | 0 | 0 | |
Values of coefficients using coefficients :
[A] = [ 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] = [ 0; 0; 1; 0; 0; 0];
X = rats(inv(A)*B)
A =
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
0 1.0000 2.0000 3.0000 4.0000 5.0000
0 0.5000 2.0000 4.5000 8.0000 12.5000
0 0.1667 1.3333 4.5000 10.6667 20.8333
0 0.0417 0.6667 3.3750 10.6667 26.0417
0 0.0083 0.2667 2.0250 8.5333 26.0417
X =
' 15/4 '
' -77/6 '
' 107/6 '
' -13 '
' 61/12 '
' -5/6 '
>> Final Equation : `(45f(i)-154f(i+1)+214f(i+2)-156f(i+3)+61f(i+4)-10f(i+5))/(12dx^2`
Central difference
The order of accuracy P = Number of nodes in the stencil N - order of derivative Q + 1
N = P+Q-1 = 5
Central difference table:
| f(i) | `Deltaxf'(i)` | `Deltax^2f^2(i)` | `Deltax^3f^3(i)` | `Deltax^4f^4(i)` | |
| af(i-2) | a | -2a | 4a/2 | -8a/6 | 16a/24 |
| bf(i-1) | b | -b | b/2 | -b/6 | b/24 |
| cf(i) | c | 0 | 0 | 0 | 0 |
| df(i+1) | d | d | d/2 | d/6 | d/24 |
| ef(i+2) | e | 2e | 4e/2 | 8e/6 | 16e/24 |
| Sum | 0 | 0 | 1 | 0 | 0 |
Values from Matlab
P = [ 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];
Q = [0;0;1;0;0];
Y = rats(inv(P)*Q)
Y =
' -1/12 '
' 4/3 '
' -5/2 '
' 4/3 '
' -1/12 'final equation =`(-f(i-2)+16f(i-1)-30f(i)+16f(i+1)-f(i+2))/(12dx^2`
Matlab main program :
clear all
close all
clc
% Given conditions
x = pi/3;
dx = linspace(pi/4, pi/400, 100);
% Analytical Function = exp(x)*cos(x)
Analytical_derivative= -2*exp(x)*sin(x);
% Taylor series Function at various node points
for i = 1:length(dx)
a(i) = (exp(x-(5*dx(i)))*cos(x-(5*dx(i))));
b(i) = (exp(x-(4*dx(i)))*cos(x-(4*dx(i))));
c(i) = (exp(x-(3*dx(i)))*cos(x-(3*dx(i))));
d(i) = (exp(x-(2*dx(i)))*cos(x-(2*dx(i))));
e(i) = (exp(x-dx(i))*cos(x-dx(i)));
f(i) = exp(x)*cos(x);
g(i) = (exp(x+dx(i))*cos(x+dx(i)));
h(i) = (exp(x+(2*dx(i)))*cos(x+(2*dx(i))));
j(i) = (exp(x+(3*dx(i)))*cos(x+(3*dx(i))));
k(i) = (exp(x+(4*dx(i)))*cos(x+(4*dx(i))));
l(i) = (exp(x+(5*dx(i)))*cos(x+(5*dx(i))));
%Central Difference Scheme
central_difference_scheme(i) = (-1*(d(i))+16*(e(i))-30*(f(i))+16*(g(i))-(h(i)))/(12*(dx(i)^2));
% Skewed left difference scheme
skewed_left_difference(i) = (45*(f(i))-154*(e(i))+214*(d(i))-156*(c(i))+61*(b(i))-10*(a(i)))/(12*(dx(i)^2));
% Skewed right difference scheme
skewed_right_difference(i) = (45*(f(i))-154*(g(i))+214*(h(i))-156*(j(i))+61*(k(i))-10*(l(i)))/(12*(dx(i)^2));
% Finding Absolute Error
central_error(i) = abs(central_difference_scheme(i)- Analytical_derivative);
skewed_left_error(i) = abs(skewed_left_difference(i)- Analytical_derivative);
skewed_right_error(i) = abs(skewed_right_difference(i)- Analytical_derivative);
end
% Plotting error vs dx values
figure(1)
loglog(dx, central_error, 'r', 'linewidth', 1.5, 'marker', 'square', 'markersize',3)
hold on
loglog(dx, skewed_left_error, 'g', 'linewidth', 1.5, 'marker', 'square', 'markersize',3)
hold on
loglog(dx, skewed_right_error , 'b', 'linewidth', 1.5, 'marker', 'square', 'markersize',3)
xlabel('dx')
ylabel('Absolute Error')
title('Comparison between Central Differencing, Skewed left differencing and Skewed right differencing schemes error')
legend('Central Differencing', 'Skewed left differencing', 'Skewed right differencing', 'location', 'southeast')
hold off
Output :
>> central_error(i)
ans =
8.4233e-10
>> skewed_left_error(i)
ans =
5.5664e-08
>> skewed_right_error(i)
ans =
5.8609e-08
From all of the observation, it is concluded that the value of the central difference scheme is very less as compared to forward as well as backward difference error.
Skewed schemes for a derivative at a point can be derived with higher orders of accuracy using additional grid points to one side of that point.