Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
AIM:
Derive the following 4th order approximations of the second-order derivative.
1. Central difference
2. Skewed right-sided difference
3. Skewed left-sided difference
After deriving these schemes, prove that your skewed schemes are fourth-order accurate
Writing a program in Matlab to evaluate the second-order derivative of the analytical function `e^(x)*cos(x)` and compare it with the 3 numerical approximations that you have derived.
Provide a plot that compares the absolute error between the above-mentioned schemes
THEORY:
TAYLOR TABLE FOR CENTRAL DIFFERENCE:
Number of Nodes = p+q-1
q = Order of Derivative = 2
p = Order of Approximation = 4
Number of Nodes/points = 2+4-1=5
`trianglex^2(d^2u)/dx^2 = af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`
| f(i) | `trianglexf'(i)` | `trianglex^2f''(x)` | `trianglex^3f'''(x)` | `trianglex^4f''''(x)` | |
| 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` |
| 0 | 0 | 1 | 0 | 0 |
From the above Taylor Table, the following algebraic equation can be formed:
a+b+c+d+e = 0
-2a-b+d+2e=0
`(4a)/2+b/2+d/2+(4e)/2 = 1`
`(-8a)/6+(-b)/6+d/6+(8e)/6 = 0`
`(16a)/24+b/24+d/24+(16e)/24 = 0`
This can be solved using the Matrix form AX=B, X= `A^-1 B`
TAYLOR TABLE FOR SKEWED RIGHT-SIDED DIFFERENCE:
Number of Nodes = p+q
q = Order of Derivative = 2
p = Order of Approximation = 4
Number of Nodes/points = 2+4=6
`trianglex^2(d^2u)/dx^2 = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)`
| f(i) | `trianglexf'(i)` | `trianglex^2f''(x)` | `trianglex^3f'''(x)` | `trianglex^4f''''(x)` | `trianglex^5f'''''(x)` | |
| 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` |
| 0 | 0 | 1 | 0 | 0 | 0 |
From the above Taylor Table, the following algebraic equation can be formed:
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`
This can be solved using the Matrix form AX=B, X= `A^-1 B`
TAYLOR TABLE FOR SKEWED LEFT-SIDED DIFFERENCE:
Number of Nodes = p+q
q = Order of Derivative = 2
p = Order of Approximation = 4
Number of Nodes/points = 2+4 = 6
`trianglex^2(d^2u)/dx^2 = af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)`
| f(i) | `trianglexf'(i)` | `trianglex^2f''(x)` | `trianglex^3f'''(x)` | `trianglex^4f''''(x)` | `trianglex^5f'''''(x)` | |
| af(i-5) | a | -5a | `(25a)/2` | `(-125a)/6` | `(625a)/24` | `(-3125a)/120` |
| bf(i-4) | b | -4b | `(16b)/2` | `(-64b)/6` | `(256b)/24` | `(-1024b)/120` |
| cf(i-3) | c | -3c | `(9c)/2` | `(-27c)/6` | `(81c)/24` | `(-243c)/120` |
| df(i-2) | d | -2d | `(4d)/2` | `(-8d)/6` | `(16d)/24` | `(-32d)/120` |
| ef(i-1) | e | -e | `e/2` | `-e/6` | `e/24` | `-e/120` |
| gf(i) | g | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 1 | 0 | 0 | 0 |
From the above Taylor Table, the following algebraic equation can be formed:
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`
This can be solved using the Matrix form AX=B, X= `A^-1 B`
FUNCTION:
`f(x) = e^x*cos(x)`
`f''(x) = -2*e^x*sin(x)`
x= `pi/4`
`trianglex =` linspace(`pi/4,pi/40,40`)
MATLAB CODE:
clear all
close all
clc
%finding the coefficients to the 4th order approximations
%central difference method
A = [1,1,1,1,1;-2,-1,0,1,2;4,1,0,1,4;-8,-1,0,1,8;16,1,0,1,16];
P = [0;0;2;0;0];
X1 = inv(A)*P
%skewed right side diffrence method
B = [1,1,1,1,1,1;0,1,2,3,4,5;0,1,4,9,16,25;0,1,8,27,64,125;0,1,16,81,256,625;0,1,32,243,1024,3125];
Q = [0;0;2;0;0;0];
X2=inv(B)*Q
%skewed LEFT side diffrence method
C = [1,1,1,1,1,1;5,4,3,2,1,0;25,16,9,4,1,0;125,64,27,8,1,0;625,256,81,16,1,0;3125,1024,243,32,1,0];
X3 = inv(C)*Q
% function input
x=pi/4;
dx=linspace(pi/4,pi/40,40);
y = exp(x)*cos(x);
exact_derivative = -2*exp(x)*sin(x)
% derivative w.r.t cental diffrence method
for i = 1:length(dx)
DR1(i) = ((X1(1)*exp(x-(2*dx(i)))*cos(x-(2*dx(i))))+(X1(2)*exp(x-dx(i))*cos(x-dx(i)))+(X1(3)*exp(x)*cos(x))+(X1(4)*exp(x+dx(i))*cos(x+dx(i)))+(X1(5)*exp(x+(2*dx(i)))*cos(x+(2*dx(i)))))/dx(i)^2;
central_error(i)=abs(exact_derivative-DR1(i));
end
% derivative w.r.t skewed right side diffrence method
for i = 1:length(dx)
DR2(i)=((X2(1)*exp(x)*cos(x))+(X2(2)*exp(x+dx(i))*cos(x+dx(i)))+(X2(3)*exp(x+(2*dx(i)))*cos(x+(2*dx(i))))+(X2(4)*exp(x+(3*dx(i)))*cos(x+(3*dx(i))))+(X2(5)*exp(x+(4*dx(i)))*cos(x+(4*dx(i))))+(X2(6)*exp(x+(5*dx(i)))*cos(x+(5*dx(i)))))/dx(i)^2;
Skewed_right_error(i) = abs(exact_derivative-DR2(i));
end
% derivative w.r.t skewed left side diffrence method
for i = 1:length(dx)
DR3(i)= ((X3(1)*exp(x-(5*dx(i)))*cos(x-(5*dx(i))))+(X3(2)*exp(x-(4*dx(i)))*cos(x-(4*dx(i))))+(X3(3)*exp(x-(3*dx(i)))*cos(x-(3*dx(i))))+(X3(4)*exp(x-(2*dx(i)))*cos(x-(2*dx(i))))+(X3(5)*exp(x-dx(i))*cos(x-dx(i)))+(X3(6)*exp(x)*cos(x)))/dx(i)^2;
Skewed_left_error(i) = abs(exact_derivative-DR3(i));
end
%comparison plots
loglog(dx,central_error,dx,Skewed_right_error,dx,Skewed_left_error)
xlabel('dx')
ylabel('Truncation Errors')
legend('Central Difference Error','Skewed Right Sided Error','Skewed Left Sided Error')
title('Variation of Absolute Error of Diffrent Schemes With Respect To The Variation In dx')
OBSERVATION/OUTPUT:
CONCLUSION:
- Variation of error with respect to decrease in dx:
Central Difference < Skewed Left Side ~< Skewed Right Side
- We can observe that central defference method produces the least error while applying fourth order accuracy.
- While considering application of CFD, in the fluid boundary region one cannot use central diffrence to compute the values as that may require points or nodes outside the domain which may still be unknown. In those cases while while computing values at mesh boundaries or regions skewed scheme is much more useful compared to central diffrence.