MATLAB - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Problem Statement: -
For the Function `f(x) = e^x*cos(x)` derive the 4th Order Accurate Approximation of 2nd Order Derivative for the given function using Taylor Table Method for the following difference Methods: -
1. Central Difference 2. Skewed right-sided difference 3. Skewed left-sided difference
Write A MATLAB Code to Solve the function using the derived 4th Order Approximate 2nd Order Derivative and compare it with the 2nd Order Derivative of Analytical Solution
Solution: -
1. AIM: -
A. To Derive 4th Order Approximation for 2nd Order Derivative using: -
a. Central Difference Method
b. Skewed Right-Side Difference Method
c. Skewed Left-Side Difference Method
B. Prove each Method to Be 4th Order Accurate
C. Write a MATLAB Code to Solve a function using all 3 Methods
D. Comment on the usefulness of Skewed Methods
2. Given: -
A. P = N - Q + ({1 for CDM} or {0 for Skewed Method})
Here,
P - Order of Accuracy of Approximation
N - Number of Nodes in Numerical Stensil
Q - Order of Derivative
CDM - Central Difference Method
B. `" Function" rarr f(x) = e^x*cos(x)`
3. Assumed: -
A. `"Differentiation Point"(x) = pi/3`
4. Analytical Derivative Formula: -
`bb ("Second Derivative" rarr (del^2f(x))/(delx^2) =- 2e^x*sin(x))`
5. Formulas & Derivation of The Approximation Methods: -
A. Central Difference Method: -
For calculating the 4th Order Approximation and 2nd Order Derivative
(P=4) = N - (Q=2) + 1
=> Number of Nodes in Numerical Stensil(N) = 5
So,
Numerical Stensil for Central Difference Method
`"Let, " (dx)^2(del^2f(x))/(delx^2) = af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`
Using Taylor Table Method For Central Difference Method
| `"Numerical Stensil " darr` | `f(i)` | `(dx)(delf(i))/(delx)` | `(dx)^2(del^2f(i))/(delx^2)` | `(dx)^3(del^3f(i))/(delx^3)` | `(dx)^4(del^4f(i))/(delx^4)` |
| `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 = (dx)^2(del^2f(x))/(delx^2)` | 0 | 0 | 1 | 0 | 0 |
On Solving,
`a = e = -1/12`
`b = d = 4/3`
`c = 5/2`
Hence,
`(dx)^2(del^2f(x))/(delx^2) = -1/12f(i-2)+4/3f(i-1)-5/2f(i)+4/3f(i+1)-1/12f(i+2)`
4th Order Approximation for 2nd Order Derivative of Central Difference Method
| `(del^2f(x))/(delx^2) = (-1/12f(i-2)+4/3f(i-1)-5/2f(i)+4/3f(i+1)-1/12f(i+2))/(dx)^2` |
B. Skewed Right-Side Difference Method: -
For Calculating 4th Order Approximation and 2nd Order Derivative
(P=4) = N - (Q=2) + 0
=> Number of Nodes in Numerical Stensil(N) = 6
So,
Numerical Stensil for Skewed Right-Side Difference Method
`"Let, " (dx)^2(del^2f(x))/(delx^2) = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+"f"f(i+5)`
| `"Numerical Stensil " darr` | `f(i)` | `(dx)(delf(i))/(delx)` | `(dx)^2(del^2f(i))/(delx^2)` | `(dx)^3(del^3f(i))/(delx^3)` | `(dx)^4(del^4f(i))/(delx^4)` | `(dx)^5(del^5f(i))/(delx^5)` |
| `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` | `(729d)/120` |
| `ef(i+4)` | `e` | `4e` | `(16e)/2` | `(64e)/6` | `(256e)/24` | `(1024e)/120` |
| `"f"f(i+5)` | `"f"` | `5"f"` | `(25"f")/2` | `(125"f")/6` | `(625"f")/24` | `(3125"f")/120` |
| `Sum = (dx)^2(del^2f(x))/(delx^2)` | 0 | 0 | 1 | 0 | 0 | 0 |
On Solving,
`a = 15/4`
`b = -77/6`
`c = 107/6`
`d = -13`
`e = 671/132`
`"f" = -5/6`
Hence,
`(dx)^2(del^2f(x))/(delx^2) = 15/4f(i)-77/6f(i+1)+107/6f(i+2)-13f(i+3)+671/132f(i+4)-5/6f(i+5)`
4th Order Approximation for 2nd Order Derivative of Skewed Right-Side Difference Method
| `(del^2f(x))/(delx^2) = (15/4f(i)-77/6f(i+1)+107/6f(i+2)-13f(i+3)+671/132f(i+4)-5/6f(i+5))/(dx)^2` |
C. Skewed Left-Side Difference Method: -
For Calculating 4th Order Approximation and 2nd Order Derivative
(P=4) = N - (Q=2) + 0
=> Number of Nodes in Numerical Stensil(N) = 6
So,
Numerical Stensil for Skewed Right-Side Difference Method
`"Let, " (dx)^2(del^2f(x))/(delx^2) = af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+"f"f(i)`
| `"Numerical Stensil " darr` | `f(i)` | `(dx)(delf(i))/(delx)` | `(dx)^2(del^2f(i))/(delx^2)` | `(dx)^3(del^3f(i))/(delx^3)` | `(dx)^4(del^4f(i))/(delx^4)` | `(dx)^5(del^5f(i))/(delx^5)` |
| `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` | `-(729c)/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` |
| `"f"f(i)` | `"f"` | 0 | 0 | 0 | 0 | 0 |
| `Sum = (dx)^2(del^2f(x))/(delx^2)` | 0 | 0 | 1 | 0 | 0 | 0 |
On Solving,
`a = -5/6`
`b = 671/132`
`c = -13`
`d = 107/6`
`e = -77/6`
`"f" = 15/4`
Hence,
`(dx)^2(del^2f(x))/(delx^2) = -5/6f(i-5)+671/132f(i-4)-13f(i-3)+107/6f(i-2)-77/6f(i-1)+15/4f(i)`
4th Order Approximation for 2nd Order Derivative of Skewed Right-Side Difference Method
| `(del^2f(x))/(delx^2) = (-5/6f(i-5)+671/132f(i-4)-13f(i-3)+107/6f(i-2)-77/6f(i-1)+15/4f(i))/(dx)^2` |
6. Proving That the Above Derived Formulas are 4th Order Accurate: -
A. Central Difference Method: -
As form proof the Coefficients of `f(x) , (dx)(delf(x))/(delx) , (dx)^3(del^3f(x))/(del^3x) and (dx)^4(del^4f(x))/(delx^4)` are equated to Zero(0) now for the Given Value of a, b, c, d, and e
Checking the Coefficient of`(dx)^5(del^5f(x))/(delx^5)`
Coefficient of `(dx)^5(del^5f(x))/(delx^5) = -(32a)/120-(b)/120+0+(d)/120+(32e)/120`
`"For", a = e = -1/12 and b = d = 4/3`
Coefficient of `(dx)^5(del^5f(x))/(delx^5) = 0`
Hence,
`.... + ("coefficient of "f^(IV)(x))(dx)^6(del^6f(x))/(delx^6) + (1)(dx)^2(del^2f(x))/(delx^2) = -1/12f(i-2)+4/3f(i-1)-5/2f(i)+4/3f(i+1)-1/12f(i+2)`
Dividing Both-sides from`(dx)^2`
`=> .... + ("coefficient of "f^(IV)(x))(dx)^4(del^6f(x))/(delx^6) + (del^2f(x))/(delx^2) = (-1/12f(i-2)+4/3f(i-1)-5/2f(i)+4/3f(i+1)-1/12f(i+2))/(dx)^2`
Ignoring All multiples of `bb (dx)^4` (=> The Order of Accuracy is 4th Order)
`(del^2f(x))/(delx^2) = (-1/12f(i-2)+4/3f(i-1)-5/2f(i)+4/3f(i+1)-1/12f(i+2))/(dx)^2`
B. Skewed Right-Side Difference Method: -
As form proof the Coefficients of `f(x) , (dx)(delf(x))/(delx) , (dx)^3(del^3f(x))/(del^3x) , (dx)^4(del^4f(x))/(delx^4) and (dx)^5(del^5f(x))/(delx^5)` are equated to Zero(0)
Hence,
`.... + ("coefficient of "f^(IV)(x))(dx)^6(del^6f(x))/(delx^6) + (1)(dx)^2(del^2f(x))/(delx^2) = 15/4f(i)-77/6f(i+1)+107/6f(i+2)-13f(i+3)+671/132f(i+4)-5/6f(i+5)`
Dividing Both-sides from`(dx)^2`
`=> .... + ("coefficient of "f^(IV)(x))(dx)^4(del^6f(x))/(delx^6) + (del^2f(x))/(delx^2) = (15/4f(i)-77/6f(i+1)+107/6f(i+2)-13f(i+3)+671/132f(i+4)-5/6f(i+5))/(dx)^2`
Ignoring All multiples of `bb (dx)^4` (=> The Order of Accuracy is 4th Order)
`(del^2f(x))/(delx^2) = (15/4f(i)-77/6f(i+1)+107/6f(i+2)-13f(i+3)+671/132f(i+4)-5/6f(i+5))/(dx)^2`
C. Skewed Left-Side Difference Method: -
As form proof the Coefficients of `f(x) , (dx)(delf(x))/(delx) , (dx)^3(del^3f(x))/(del^3x) , (dx)^4(del^4f(x))/(delx^4) and (dx)^5(del^5f(x))/(delx^5)` are equated to Zero(0)
Hence,
`.... + ("coefficient of "f^(IV)(x))(dx)^6(del^6f(x))/(delx^6) + (1)(dx)^2(del^2f(x))/(delx^2) = -5/6f(i-5)+671/132f(i-4)-13f(i-3)+107/6f(i-2)-77/6f(i-1)+15/4f(i)`
Dividing Both-sides from`(dx)^2`
`=> .... + ("coefficient of "f^(IV)(x))(dx)^4(del^6f(x))/(delx^6) + (del^2f(x))/(delx^2) = (-5/6f(i-5)+671/132f(i-4)-13f(i-3)+107/6f(i-2)-77/6f(i-1)+15/4f(i))/(dx)^2`
Ignoring All multiples of `bb (dx)^4` (=> The Order of Accuracy is 4th Order)
`(del^2f(x))/(delx^2) = (-5/6f(i-5)+671/132f(i-4)-13f(i-3)+107/6f(i-2)-77/6f(i-1)+15/4f(i))/(dx)^2`
7. Explanation of Code: -
A. Section A(Defining Variables): -
a. Defined x and a vector of different values of dx
b. Calculated Analytical Value using Analytical Formula
c. created "A" and "B" of Central and Skewed Forms ( "A" and "B" of AX = B)
d. Calculated "X" of Central and Skewed Forms using X = AB
B. Section B(Differentiating): -
a. Created A loop of Iteration of the length of dx vector
b. Selected a dx at each iteration
c. Created small functions for simplification of Calculations at each iteration
d. Calculated Central Difference and Skewed Differences using the Function defined above and X calculated in Section A.
e. Calculated Error generated per Iteration w.r.t Analytical Value Calculated in Section A and saved it in a Vector.
C. Section C(Plotting)
a. Plotted a log Graph between The Value of Error vs The Value of dx for all 3 Methods
b. Defined Proper Legends and Labels in the plot.
8. Code: -
% 4th Order 2nd Differentiation using Central and Both Skewed Schemes
clear all
close all
clc
%% Defining Variables
x = pi/3;
dx_all = x.*(10.^[-3:-1]);
analytical_val = -(2*exp(x)*sin(x))
A_crt_vec = [-2:1:2]';
A_rgt_vec = [0:5]';
A_lft_vec = [-5:0]';
A_crt = [ones(5,1),(A_crt_vec.^[1:4]).*(1./cumprod([1:4]))]';
A_rgt = [ones(6,1),(A_rgt_vec.^[1:5]).*(1./cumprod([1:5]))]';
A_lft = [ones(6,1),(A_lft_vec.^[1:5]).*(1./cumprod([1:5]))]';
B_crt = [0;0;1;0;0];
B_skw = [0;0;1;0;0;0];
X_crt = A_crtB_crt
X_rgt = A_rgtB_skw
X_lft = A_lftB_skw
%% Differentiating
for i=1:length(dx_all)
% Selecting h and Creating small Function for Ease
dx = dx_all(i);
f_mi_5 = exp(x-5*dx)*cos(x-5*dx);
f_mi_4 = exp(x-4*dx)*cos(x-4*dx);
f_mi_3 = exp(x-3*dx)*cos(x-3*dx);
f_mi_2 = exp(x-2*dx)*cos(x-2*dx);
f_mi_1 = exp(x-dx)*cos(x-dx);
f_x = exp(x)*cos(x);
f_pl_1 = exp(x+dx)*cos(x+dx);
f_pl_2 = exp(x+2*dx)*cos(x+2*dx);
f_pl_3 = exp(x+3*dx)*cos(x+3*dx);
f_pl_4 = exp(x+4*dx)*cos(x+4*dx);
f_pl_5 = exp(x+5*dx)*cos(x+5*dx);
% Creating Formula and Executing it
crt_fn_vec = [f_mi_2,f_mi_1,f_x,f_pl_1,f_pl_2];
skw_rgt_fn_vec = [f_x,f_pl_1,f_pl_2,f_pl_3,f_pl_4,f_pl_5];
skw_lft_fn_vec = [f_mi_5,f_mi_4,f_mi_3,f_mi_2,f_mi_1,f_x];
dif_crt = (crt_fn_vec*X_crt)/(dx^2);
dif_rgt = (skw_rgt_fn_vec*X_rgt)/(dx^2);
dif_lft = (skw_lft_fn_vec*X_lft)/(dx^2);
% Storing Errors
err_crt(i) = abs(analytical_val - dif_crt);
err_rgt(i) = abs(analytical_val - dif_rgt);
err_lft(i) = abs(analytical_val - dif_lft);
end
%% Plotting OUTPUT
loglog(dx_all,err_crt,dx_all,err_rgt,dx_all,err_lft)
xlabel('Value of dx','FontSize',15)
ylabel('Value of Error','FontSize',15)
legend('Central Difference','Skewed Right-Side Difference','Skewed Left-Side Difference')9. OUTPUTS Generated on Running the Code: -
A. Command Window
B. Plot OUTPUT
10. Usefulness of Skewed Difference Methods: -
Skewed Schemes are Generally useful where the Information is flowing from one direction then the Correct calculations can only be developed by points on one side, as the next value thus defined is required to be developed using the values just achieved. To explain this I'm using the example Below.
Example,
In this example, if properties at points A and B are available and the Value of C is to be determined such that difference between each point is dx then,
Using Taylor Series, `phi_C = phi_B + ((delphi)/(delx))_Bdx`
Now, As properties at points A and B are only Known So to determine `((delphi)/(delx))_B`
Skewed Left-Side Difference Method is used, `((delphi)/(delx))_B = (phi_B - phi_A)/dx`