Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Derivation:
Consider the second derivative of function f(x) be d2f/dx2.
- Using Central Differencing Scheme
First we need to determine the order of accuracy needed:
The order of accuracy P = Number of nodes in the stencil N - order of derivative Q + 1
Hence N = P + Q - 1
N = 2+4-1
N=5
Considering the 5 nodes as shown below for CD scheme
Then, d2f/dx2 = (af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2
Applying taylor series expansion as follow:
a*f(i-2) = a[f(i)-f’(i)(2dx)+(f’’(i)*(2dx)^2)/2!-((f’’’(i)*(2dx)^3)/3!)+((f’’’’(i)*(2dx)^4)/4!]
b*f(i-1) = b[f(i)-f’(i)dx+(f’’(i)*dx^2)/2!-((f’’’(i)*dx^3)/3!)+((f’’’’(i)*dx^4)/4!]
c*f(i) = c*f(i)
d*f(i+1) = d[f(i)+f’(i)dx+(f’’(i)*dx^2)/2!+((f’’’(i)*dx^3)/3!)+((f’’’’(i)*dx^4)/4!]
e*f(i+2) = e[f(i)+f’(i)(2dx)+(f’’(i)*(2dx)^2)/2!+((f’’’(i)*(2dx)^3)/3!)+((f’’’’(i)*(2dx)^4)/4!]
|
f(i) |
f’(i)*dx |
f’’(i)*dx^2 |
f’’’(i)*dx^3 |
f’’’’(i)*dx^4 |
|
|
a*f(i-2) |
a |
-2a |
2a |
-4a/3 |
2a/3 |
|
b*f(i-1) |
b |
-b |
b/2 |
-b/6 |
b/24 |
|
c*f(i) |
c |
0 |
0 |
0 |
0 |
|
d*f(i+1) |
d |
d |
d/2 |
d/6 |
d/24 |
|
e*f(i+2) |
e |
2e |
2e |
4e/3 |
2e/3 |
|
Taylor table method |
0 |
0 |
1 |
0 |
0 |
Hence our system of equations will be
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
The same in the matrix form will be
A = [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]
X = [a; b; c; d; e]
B = [0; 0; 1; 0; 0]
Solving these equations,
a = -1/12
b = 4/3
c = -5/2
d = 4/3
e = -1/12
Applying these constant values back to second order equation,
d2f/dx2 = (af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2
d2f/dx2 = (-1/12)*f(i-2)+(4/3)*f(i-1)-(5/2)*f(i)+(4/3)*f(i+1)-(1/12)*f(i+2)/dx^2
- Using Skewed Right Sided Differencing Scheme
First we need to determine the order of accuracy needed:
The order of accuracy P = Number of nodes in the stencil N - order of derivative Q
Hence N = P + Q
N = 2+4
N=6
Considering the 6 nodes as shown below for Skewed RSD scheme
Then, d2f/dx2 = (af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/dx^2
Applying taylor series expansion as follow:
af(i) = a*f(i)
b*f(i+1) = b[f(i)+f’(i)dx+(f’’(i)*dx^2)/2!+((f’’’(i)*dx^3)/3!)+((f’’’’(i)*dx^4)/4!+((f’’’’’(i)*dx^5)/5!)]
c*f(i+2) = c[f(i)+f’(i)(2dx)+(f’’(i)*(2dx)^2)/2!+((f’’’(i)*(2dx)^3)/3!)+((f’’’’(i)*(2dx)^4)/4!)+((f’’’’’(i)*(2*dx)^5)/5!)]
d*f(i+3) = d[f(i)+f’(i)(3dx)+(f’’(i)*(3dx)^2)/2!+((f’’’(i)*(3dx)^3)/3!)+((f’’’’(i)*(3dx)^4)/4!)+((f’’’’’(i)*(3*dx)^5)/5!)]
e*f(i+4) = e[f(i)+f’(i)(4dx)+(f’’(i)*(4dx)^2)/2!+((f’’’(i)*(4dx)^3)/3!)+((f’’’’(i)*(4dx)^4)/4!)+((f’’’’’(i)*(4*dx)^5)/5!)]
g*f(i+5) = g[f(i)+f’(i)(5dx)+(f’’(i)*(5dx)^2)/2!+((f’’’(i)*(5dx)^3)/3!)+((f’’’’(i)*(5dx)^4)/4!)+((f’’’’’(i)*(5*dx)^5)/5!)]
|
f(i) |
f’(i)*dx |
f’’(i)*dx^2 |
f’’’(i)*dx^3 |
f’’’’(i)*dx^4 |
f’’’’’(i)*dx^5 |
|
|
a*f(i) |
a |
0 |
0 |
0 |
0 |
0 |
|
b*f(i+1) |
b |
b |
b/2 |
b/6 |
b/24 |
b/120 |
|
c*f(i+2) |
c |
2c |
4c/2 |
8c/6 |
16c/24 |
32c/120 |
|
d*f(i+3) |
d |
3d |
9d/2 |
27d/6 |
81d/24 |
243d/120 |
|
e*f(i+4) |
e |
4e |
16e/2 |
64e/6 |
256e/24 |
1024e/120 |
|
g*f(i+5) |
g |
5g |
25g/2 |
125g/6 |
625g/24 |
3125g/120 |
|
Taylor table method |
0 |
0 |
1 |
0 |
0 |
0 |
Hence our system of equations will be
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 same in the matrix form will be
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]
X = [a; b; c; d; e; g]
B = [0; 0; 1; 0; 0; 0]
Solving these equations,
a = 15/4
b = -77/6
c = 107/6
d = -13
e = 61/12
f = -5/6
Applying these constant values back to second order equation,
d2f/dx2 = (af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/dx^2
d2f/dx2 = (15/4)*f(i)-(77/6)*f(i+1)+(107/6)*f(i+2)-13*f(i+3)+(61/12)*f(i+4)-(5/6)*f(i+5)/dx^2
- Using Skewed Left Sided Differencing Scheme
First we need to determine the order of accuracy needed:
The order of accuracy P = Number of nodes in the stencil N - order of derivative Q
Hence N = P + Q
N = 2+4
N=6
Considering the 6 nodes as shown below for Skewed LSD scheme
Then, d2f/dx2 = (af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i))/dx^2
Applying taylor series expansion as follow:
a*f(i-5) = a*[f(i)-f’(i)(5dx)+(f’’(i)*(5dx)^2)/2!-((f’’’(i)*(5dx)^3)/3!)+((f’’’’(i)*(5dx)^4)/4!-((f’’’’’(i)*(5dx)^5)/5!)]
b*f(i-4) = b*[f(i)-f’(i)(4dx)+(f’’(i)*(4dx)^2)/2!-((f’’’(i)*(4dx)^3)/3!)+((f’’’’(i)*(4dx)^4)/4!-((f’’’’’(i)*(4dx)^5)/5!)]
c*f(i-3) = c*[f(i)-f’(i)(3dx)+(f’’(i)*(3dx)^2)/2!-((f’’’(i)*(3dx)^3)/3!)+((f’’’’(i)*(3dx)^4)/4!-((f’’’’’(i)*(3dx)^5)/5!)]
d*f(i-2) = d*[f(i)-f’(i)(2dx)+(f’’(i)*(2dx)^2)/2!-((f’’’(i)*(2dx)^3)/3!)+((f’’’’(i)*(2dx)^4)/4!-((f’’’’’(i)*(2dx)^5)/5!)]
e*f(i-1) = e*[f(i)-f’(i)dx+(f’’(i)*dx^2)/2!-((f’’’(i)*dx^3)/3!)+((f’’’’(i)*dx^4)/4!-((f’’’’’(i)*dx^5)/5!)]
gf(i) = g*f(i)
|
f(i) |
f’(i)*dx |
f’’(i)*dx^2 |
f’’’(i)*dx^3 |
f’’’’(i)*dx^4 |
f’’’’’(i)*dx^5 |
|
|
a*f(i-5) |
a |
-5a |
25a/2 |
-125a/6 |
625a/24 |
-3125a/120 |
|
b*f(i-4) |
b |
-4b |
16b/2 |
-64b/6 |
256b/24 |
-1024b/120 |
|
c*f(i-3) |
c |
-3c |
9c/2 |
-27c/6 |
81c/24 |
-243c/120 |
|
d*f(i-2) |
d |
-2d |
4d/2 |
-8d/6 |
16d/24 |
-32d/120 |
|
e*f(i-1) |
e |
-e |
e/2 |
-e/6 |
e/24 |
-e/120 |
|
g*f(i) |
g |
0 |
0 |
0 |
0 |
0 |
|
Taylor table method |
0 |
0 |
1 |
0 |
0 |
0 |
Hence our system of equations will be
a+b+c+d+e+g=0
-5a-4b-3c-2d-e=0
(25a/2)+(16b/2)+(9c/2)+(4d/2)+(e/2)=1
(-125a/6)-(64b/6)-(27c/6)-(8d/6)-(e/6)=0
(625a/24)+(256b/24)+(81c/24)+(16d/24)+(e/24)=0
(-3125a/120)-(1024b/120)-(243c/120)-(32d/120)-(e/120)=0
The same in the matrix form will be
A = [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]
X = [a; b; c; d; e; g]
B = [0; 0; 1; 0; 0; 0]
Solving these equations,
a = -5/6
b = 61/12
c = -13
d = 107/6
e = -77/6
g = 15/4
Applying these constant values back to second order equation,
d2f/dx2 = (af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i))/dx^2
d2f/dx2 = -(5/6)*f(i-5)+(61/12)*f(i-4)-13*f(i-3)+(107/6)*f(i-2)-(77/6)*f(i-1)+(15/4)*f(i))/dx^2
==============================================================
Analytical derivative
Let f(x) = exp(x)*cos(x)
f’(x) = exp(x)*cos(x)-sin(x)*exp(x)
f’’(x) = exp(x)*cos(x)-sin(x)*exp(x)-exp(x)*cos(x)-sin(x)*exp(x)
f’’(x) = -2*exp(x)*sin(x)
Describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?
Skewed schemes used in finite difference method uses non-equidistant points to approximate the derivative of the function. This technique is very useful where boundary conditions are not compatible to central difference scheme.
Comparatively referring to the resultant graph, skewed scheme seems more accurate than central differencing scheme. This is due to number of points (which is more) used to derivate approximation of function and thus reduces error.
Also since the skewed scheme does not require values from both sides of the point of interest, it can be used to approximate derivatives at boundaries, which CD scheme cannot.
Matlab code is shown below:
clear all
close all
clc
% to evaluate the second-order derivative of the analytical function exp(x)*cos(x)
%analytical derivative f''(x) = -2*exp(x)*sin(x)
x = pi/3;
dx = linspace(pi/4000, pi/4, 40);
analy_deriv = -2*exp(x)*sin(x);
for i=1:length(dx)
% comparing with central differencing scheme
% Replaced decimals with fractions a = -0.0833; b = 1.3333; c = -2.5000; d = 1.3333; e = -0.0833;
a = -1/12;
b = 4/3;
c = -5/2;
d = 4/3;
e = -1/12;
% af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2
cds(i) = (a*exp(x-2*dx(i))*cos(x-2*dx(i))+(b*exp(x-dx(i))*cos(x-dx(i)))+(c*exp(x)*cos(x))+(d*exp(x+dx(i))*cos(x+dx(i)))+(e*exp(x+2*dx(i))*cos(x+2*dx(i))))/(dx(i)^2);
error_cds(i) = abs(cds(i)-analy_deriv);
% comparing with right skewed differencing scheme
% Replaced decimals with fractions a1 = 3.7500; b1 = -12.8333;c1 = 17.8333;d1 = -13.0000;e1 = 5.0833;g1 = -0.8333;
a1 = 15/4;
b1 = -77/6;
c1 = 107/6;
d1 = -13;
e1 = 61/12;
g1 = -5/6;
% (af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/dx^2
rsd(i) = (a1*exp(x)*cos(x)+(b1*exp(x+dx(i))*cos(x+dx(i)))+(c1*exp(x+2*dx(i))*cos(x+2*dx(i)))+(d1*exp(x+3*dx(i))*cos(x+3*dx(i)))+(e1*exp(x+4*dx(i))*cos(x+4*dx(i)))+(g1*exp(x+5*dx(i))*cos(x+5*dx(i))))/(dx(i)^2);
error_rsd(i) = abs(rsd(i)-analy_deriv);
% comparing with left skewed differencing scheme
% Replaced decimals with fractions a2 = -0.8333;b2 = 5.0833;c2 = -13.0000;d2 = 17.8333;e2 = -12.8333;g2 = 3.7500;
a2 = -5/6;
b2 = 61/12;
c2 = -13;
d2 = 107/6;
e2 = -77/6;
g2 = 15/4;
% (af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i))/dx^2
lsd(i) = (a2*exp(x-5*dx(i))*cos(x-5*dx(i))+(b2*exp(x-4*dx(i))*cos(x-4*dx(i)))+(c2*exp(x-3*dx(i))*cos(x-3*dx(i)))+(d2*exp(x-2*dx(i))*cos(x-2*dx(i)))+(e2*exp(x-dx(i))*cos(x-dx(i)))+(g2*exp(x)*cos(x)))/(dx(i)^2);
error_lsd(i) = abs(lsd(i)-analy_deriv);
end
loglog(dx,error_cds,'r');
hold on;
loglog(dx,error_lsd,'g');
hold on;
loglog(dx,error_rsd,'b');
% grid on;
xlabel('dx');
ylabel('Absolute Error');
title('Fourth Order Approximation');
legend('Central Differencing Scheme', 'Skewed Left Sided Differencing Scheme', 'Skewed Right Sided Differencing Scheme');