Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Aim: To derive 4th order approximation of a 2nd order derivative using Taylor Table method in MATLAB.
Objective:
- Derive the following 4th order approximations of the second-order derivative.
-
Central difference
-
Skewed right-sided difference
-
Skewed left-sided difference
-
Theory:
-
Central difference
Finding the no of nodes:
`rarr `Then the Number of nodes N = The order of accuracy P + Order of derivative Q - 1
`rarr ` `N = 4+2-1`
`rarr ` `N = 5`- Derivation:
`(del^2f)/(delx^2) = af(i-2) + bf(i-j) + cf(i) + df(i+1) +ef(i+2)` -
Taylor series expansion for the above terms:
`af(i-2) = af(i) - (2dxf^'(i))/1 + (4adx^2f^''(i))/2 - (4adx^3f^'''(i))/3 + (2adx^4f^'''(i))/3`
`bf(i-1) = bf(i) - (bdxf^'(i))/1 + (bdx^2f^''(i))/2 - (bdx^3f^'''(i))/3 + (bdx^4f^'''(i))/3`
`to cf(i) = cf(i)`
`df(i+1) = df(i) + (d.dxf'(i))/(1!) + (d.dx^2f''(i))/(2!) + (d.dx^3f'''(i))/(3!) + (d.dx^4f''''(i))/(4!)`
`ef(i+2) = ef(i) + (2edxf'(i))/(1!) + (4edx^2f''(i))/(2!) + (8edx^3f'''(i))/(3!) + (16dx^4f''''(i))/(4!)`-
Representing in the form of a table:
`dx^2f''(i)` `f(i)` `f'(i)dx` `f''(i)dx^2` `f''''(i)dx^3` `f'''''(i)dx^4` `af(i-2)` `a` `-2a` `2a` `(-4a)/3` `(2a)/3` `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` `2e `4e/3` `(2e)/3` Sum of coefficient `0` `0` `1` `0` `0` - `A = [[1,1,1,1,1], [-2,-1,0,0,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]]`
`to X = (inv)A*B`
`therefore [[1,1,1,1,1], [-2,-1,0,0,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]]` * `[[a],[b],[c],[d],[e]]` * `[[0],[0],[1],[0],[0]]` - We get the value of the unknown is
`a = -0.0833, b = 1.333, c = -2.500, d = 1.333, e = -0.0833`
-
- Derivation:
-
Right side difference:
The order of accuracy P = Number of nodes in the stencil N - order of derivative Q for one-sided/skewed schemes.
hear, `N= 6`
`(del^2f)/(delx^2) = af(i) + bf(i+1) + cf(i+2) + df(i+3) + ef(i+4) + ff(i+5)`-
Representing in the form of a table:
`f(x)` `dxf'(x)` `dx^2f''(x)` `dx^3f'''(x)` `dx^4f''''(x)` `dx^5f''''''(x)` `af(1)` `0` `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` `ff(i+5)` `f` `5f` `(25f)/2` `(125f)/6` `(625f)/24` `(3125f)/120` Sum of coefficient `0` `0` `1` `0` `0` `0`
- `A = [[1,1,1,1,1,1], [0,1,2,3,4,5], [0,1/2,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],[f]]`
- `B = [[0],[0],[1],[0],[0],[0]]`
`therefore [[1,1,1,1,1,1], [0,1,2,3,4,5], [0,1/2,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]] * [[a],[b],[c],[d],[e],[f]] * [[0],[0],[1],[0],[0],[0]]` - We get the value of the unknown is
`a = 3.750, b = -12.833, c = 17.833, d = -13, e = -0.0833`
-
- Skewed left sided difference:
`(del^2f)/(delx^2) = af(i-5) + bf(i-4) + cf(i-3) + df(i-2) + ef(i-1) + ff(i)`
- Representing in the form of a table:
`f(x)` `dxf'(x)` `dx^2f''(x)` `dx^3f'''(x)` `dx^4f'''''(x)` `dx^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` `ff(i) `f` `0` `0` `0` `0` `0`` `af(1-5)+bf(i-4)+
cf(i-3)+df(i-2)+
ef(i-1)+ff(i)`
`0` `0` `0` `0` `0`
- `A = [[1,1,1,1,1,1], [-5,-4,-3,-2,-1,0], [25/2,16/2,9/2,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],[f]]`
- `B = [[0],[0],[1],[0],[0],[0]]`
`therefore A = [[1,1,1,1,1,1], [-5,-4,-3,-2,-1,0], [25/2,16/2,9/2,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]] * [[a],[b],[c],[d],[e],[f]] * [[0],[0],[1],[0],[0],[0]]` - We get the value of the unknown is
`a = -0.833, b = 5.083, c = 13, d = 17.833, e = 3.75`
- Representing in the form of a table:
Matlab Code
- Function Code
-
Central Differencing:
function Error_central_differencing = fourth_order_central_differencing(x,dx); %analytic function = exp(x)*cos(x) %analytical derivative %f''(x)=-2*exp(x)sin(x) analytical_derivative=-2*exp(x)*sin(x); %numerical derivative %central differencing central_differencing=(((-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_differencing=abs(central_differencing-analytical_derivative); end -
Right sided skew difference:
function Error_right_sided_skew = fourth_order_right_sided_skew(x,dx); %analytic function = exp(x)*cos(x) %analytical derivative %f''(x)=-2*exp(x)sin(x) analytical_derivative=-2*exp(x)*sin(x); %numerical derivative %right sided skew right_sided_skew=(((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_sided_skew=abs(right_sided_skew-analytical_derivative); end -
Left sided skew difference:
function Error_left_sided_skew = fourth_order_left_sided_skew(x,dx); %analytic function = exp(x)*cos(x) %analytical derivative %f''(x)=-2*exp(x)sin(x) analytical_derivative=-2*exp(x)*sin(x); %numerical derivative %Left sided skew left_sided_skew=(((-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_sided_skew=abs(left_sided_skew-analytical_derivative); end - Main Code
clear all close all clc %Value of x x=pi/3; %Range of dx dx=linspace(pi/4,pi/4000,20); %For loop for i=1:length(dx) central_differencing(i)=fourth_order_central_differencing(x,dx(i)); right_sided_skew(i)=fourth_order_right_sided_skew(x,dx(i)); left_sided_skew(i)=fourth_order_left_sided_skew(x,dx(i)); end %Plotting figure(1) loglog(dx,central_differencing,'g', 'LineWidth',2) hold on loglog(dx,right_sided_skew,'b','LineWidth',2) hold on loglog(dx,left_sided_skew,'r','LineWidth',2) grid on xlabel('Range of dx') ylabel('Truncation error') title('Range of dx vs Error') legend('Central differencing','Skewed rigth sided difference','Skewed left sided difference','location','southeast');
-
Result:
Conclusion:
- As seen from the plot the central differencing scheme has very less error when compared to the other schemes. Therefore, the central differencing scheme is better than other schemes because it uses information from both the sides from the point of consideration to give appropriate and accurate values.
- Skewed right gives maximum error or in other words we can say that it's unstable as compared to other schemes.
- Central difference gives the lowest error or in other words we can say that its stable as compared to other schemes.
- Skewed left gives error somewhere in between skewed right difference and central difference.