Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
OBJECTIVE:
1.To derive the following 4th order approximations of the second-order derivative.
a. Central difference
b. Skewed right-sided difference
c. Skewed left-sided difference
2.To write a program in Matlab to evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations that you have derived.
3.To Provide a plot that compares the absolute error between the above-mentioned schemes
Taylor's Series:
The Taylor series of a function is an infinite sum of terms that are expressed in terms of the function's derivatives at a single point.
`f(x)=f(a)+(f'(a))/(1!)(x-a)+(f''(a))/(2!)(x-a)^2+(f'''(a))/(3!)(x-a)^3+...,`
a.Central differencing:
We are going to derive the central differencing for the fourth order approximation of the second order derivative of a function using Taylor's table method.
In order to derive the central differencing we need to calculate the number of nodes,
For Central differencing,Number of Nodes, N=P+Q-1
where,
P = Order of the derivative.
Q = Order of Approximation.
For fourth order approximation of the second order derivative , the number of nodes required to derive central differencing is,
N = 2 + 4 - 1 = 5
The hypothesis of the second order derivative is expressed as,
`nablax^2f''(x) = af(i-2) + bf(i-1) + cf(i) + df(i + 1) + ef(i + 2) ...........(1)`
Now evaluating each term of the eq(1) using taylor series,
`af(i - 2) = a[f(i) - (2nablax)/(1!)f'(i) + (2nablax)^2/(2!)f''(i) - (2nablax)^3/(3!)f'''(i) + (2nablax)^4/(4!) - .........]`
`bf(i - 1) = b[f(i) - (nablax)/(1!)f'(i) + (nablax)^2/(2!)f''(i) - (nablax)^3/(3!)f'''(i) + (nablax)^4/(4!) - .........]`
`cf(i ) = cf(i)`
`df(i + 1) = d[f(i) + (nablax)/(1!)f'(i) + (nablax)^2/(2!)f''(i) + (nablax)^3/(3!)f'''(i) + (nablax)^4/(4!) - .........]`
`ef(i + 2) = e[f(i) + (2nablax)/(1!)f'(i) + (2nablax)^2/(2!)f''(i) + (2nablax)^3/(3!)f'''(i) + (2nablax)^4/(4!) - .........]`
Generating Taylor's table to visualise the equations properly,
| `nablax^2(del^2f)/(delx^2)` | `f(i)` | `nablaxf'(i)` | `nablax^2f''(i)` | `nablax^3f'''(i)` | `nablax^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 |
| sum of coefficients | 0 | 0 | 1 | 0 | 0 |
With the help of the taylor's table we can find out the values of the coefficients by forming the matrix and find the values in MATLAB;
`[[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]][[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`
MATLAB SYNTAX FOR DERIVING THE VALUES
a_central = [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];
b_central = [0;0;1;0;0];
c_central = inv(a_central)*b_central;
%Solution
>> a_central
a_central =
1.0000 1.0000 1.0000 1.0000 1.0000
-2.0000 -1.0000 0 1.0000 2.0000
2.0000 0.5000 0 0.5000 2.0000
-1.3333 -0.1667 0 0.1667 1.3333
0.6667 0.0417 0 0.0417 0.6667
>> b_central
b_central =
0
0
1
0
0
>> c_central = inv(a_central)*b_central
c_central =
-0.0833
1.3333
-2.5000
1.3333
-0.0833On solving the above equations we will get the values for coeffiecients and these coefficients value will then be put in equation 1 and hence the second order derivative for the fourth order approximation of the function for central differencing is evaluated.
Skewed right-sided difference:
In this case we are going to follow the same as procedure as we did in central differencing.
The number of nodes required,
N=P+Q
P = Order of the derivative.
Q = Order of Approximation.
For fourth order approximation of the second order derivative , the number of nodes required to derive skewed right-sided difference is,
N = 2 + 4 = 6.
The hypothesis of the second order derivative is expressed as,
`nablax^2f''(x) = af(i) + bf(i+1) + cf(i+2) + df(i + 3) + ef(i + 4) + gf(i+5) ...........(2)`
Now evaluating each term of the eq(2) using taylor series,
`af(i) = af(i)`
`bf(i + 1) = b[f(i) + (nablax)/(1!)f'(i) + (nablax)^2/(2!)f''(i) + (nablax)^3/(3!)f'''(i) + (nablax)^4/(4!) - .........]`
`cf(i+2) = c[f(i) + (2nablax)/(1!)f'(i) + (2nablax)^2/(2!)f''(i) + (2nablax)^3/(3!)f'''(i) + (2nablax)^4/(4!) - .........]`
`df(i + 3) = d[f(i) + (3nablax)/(1!)f'(i) + (3nablax)^2/(2!)f''(i) + (3nablax)^3/(3!)f'''(i) + (3nablax)^4/(4!) - .........]`
`ef(i + 4) = e[f(i) + (4nablax)/(1!)f'(i) + (4nablax)^2/(2!)f''(i) + (4nablax)^3/(3!)f'''(i) + (4nablax)^4/(4!) - .........]`
`gf(i + 5) = g[f(i) + (5nablax)/(1!)f'(i) + (5nablax)^2/(2!)f''(i) + (5nablax)^3/(3!)f'''(i) + (5nablax)^4/(4!) - .........]`
| `nablax^2(del^2f)/(delx^2)` | `f(i)` | `nablaxf'(i)` | `nablax^2f''(i)` | `nablax^3f'''(i)` | `nablax^4f''''(x)` | `nablax^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 |
| Evaluation for 2nd order derivative | 0 | 0 | 1 | 0 | 0 | 0 |
With the help of the taylor's table we can find out the values of the coefficients by forming the matrix and find the values in MATLAB;
`[[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]][[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`
MATLAB SYNTAX FOR DERIVING THE VALUES
a_right=[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_right=[0;0;1;0;0;0]
c_right = inv(a_right)*b_right
%solution
>> a_right
a_right =
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
>> b_right
b_right =
0
0
1
0
0
0
>> c_right
c_right =
3.7500
-12.8333
17.8333
-13.0000
5.0833
-0.8333
>> On solving the above equations we will get the values for coeffiecients and these coefficients value will then be put in equation 2 and hence the second order derivative for the fourth order approximation of the function for right-sided differencing is evaluated.
Skewed left-sided difference:
In this case we are going to follow the same as procedure as we did in right - sided differencing.
The number of nodes required,
N=P+Q
P = Order of the derivative.
Q = Order of Approximation.
For fourth order approximation of the second order derivative , the number of nodes required to derive skewed right-sided difference is,
N = 2 + 4 = 6.
The hypothesis of the second order derivative is expressed as,
`nablax^2f''(x) = af(i-5) + bf(i-4) + cf(i-3) + df(i -2) + ef(i -1) + gf(i) ...........(3)`
Now evaluating each term of the eq(3) using taylor series,
`af(i-5) = a[f(i) - (5nablax)/(1!)f'(i) + (5nablax)^2/(2!)f''(i) - (5nablax)^3/(3!)f'''(i) + (5nablax)^4/(4!) - .........]`
`bf(i - 4) = b[f(i) - (4nablax)/(1!)f'(i) + (4nablax)^2/(2!)f''(i) - (4nablax)^3/(3!)f'''(i) + (4nablax)^4/(4!) - .........]`
`cf(i-3) = c[f(i) - (3nablax)/(1!)f'(i) + (3nablax)^2/(2!)f''(i) - (3nablax)^3/(3!)f'''(i) + (3nablax)^4/(4!) - .........]`
`df(i -2) = d[f(i) - (2nablax)/(1!)f'(i) + (2nablax)^2/(2!)f''(i) - (2nablax)^3/(3!)f'''(i) + (2nablax)^4/(4!) - .........]`
`ef(i - 1) = e[f(i) - (nablax)/(1!)f'(i) + (nablax)^2/(2!)f''(i) - (nablax)^3/(3!)f'''(i) + (nablax)^4/(4!) - .........]`
`gf(i) = gf(i)`
| `nablax^2(del^2f)/(delx^2)` | `f(i)` | `nablaxf'(i)` | `nablax^2f''(i)` | `nablax^3f'''(i)` | `nablax^4f''''(x)` | `nablax^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 |
| `bf(i-3)` | c | -3c | 9c/2 | -27c/6 | 81c/24 | -243c/120 |
| `bf(i-2)` | d | -2d | 4d/2 | -8d/6 | 16d/24 | -32d/120 |
| `bf(i-1)` | e | -e | e/2 | -e/6 | e/24 | -e/120 |
| `gf(i)` | g | 0 | 0 | 0 | 0 | 0 |
| Evaluation for 2nd order derivative | 0 | 0 | 1 | 0 | 0 | 0 |
With the help of the taylor's table we can find out the values of the coefficients by forming the matrix and find the values in MATLAB;
`[[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]][[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`
MATLAB SYNTAX FOR DERIVING THE VALUES
a_left = [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];
b_left = [0;0;1;0;0;0];
c_left = inv(a_left) * b_left;
%solution
>> a_left
a_left =
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
-5.0000 -4.0000 -3.0000 -2.0000 -1.0000 0
12.5000 8.0000 4.5000 2.0000 0.5000 0
-20.8333 -10.6667 -4.5000 -1.3333 -0.1667 0
26.0417 10.6667 3.3750 0.6667 0.0417 0
-26.0417 -8.5333 -2.0250 -0.2667 -0.0083 0
>> b_left
b_left =
0
0
1
0
0
0
>> c_left
c_left =
-0.8333
5.0833
-13.0000
17.8333
-12.8333
3.7500
>> On solving the above equations we will get the values for coeffiecients and these coefficients value will then be put in equation 3 and hence the second order derivative for the fourth order approximation of the function for left-sided differencing is evaluated
MAIN SYNTAX FOR DETERMINING THE VARIATION OF DIFFERENT SCHEME USING MATLAB
clear all
close ALL
clc
%Using taylor table method we are going to find out the second order
%derivative for fourth order approximation
% Function, f(x) = exp(x).cos(x)
x=pi/3;
dx=linspace(pi/4000,pi/40,100);
% Calculating the second order analytical derivative of the function
analytical_derivative = -2*exp(x)*sin(x);
%Central Differencing
a_central = [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];
b_central = [0; 0; 1; 0; 0];
c_central = inv(a_central)*b_central;
%Calculating second order numerical derivative using Central differencing
central_differencing = (c_central(1)*(exp(x-dx*2).*sin(x-dx*2)) + ...
c_central(2)*(exp(x-dx).*sin(x-dx)) + c_central(3)*(exp(x).*sin(x)) + ...
c_central(4).*(exp(x+dx).*sin(x+dx)) + c_central(5).*(exp(x+dx*2).*sin(x+dx*2)))./(dx.^2);
% Absolute error between second order derivative of analytical and central
% differencing scheme
central_error = abs(central_differencing - analytical_derivative);
%Right Sided Skewed
a_right = [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_right = [0;0;1;0;0;0];
c_right = inv(a_right) * b_right;
% Calculating second order numerical derivative using Central differencing
skewed_right = (c_right(1).*(exp(x)*sin(x)) + c_right(2).*(exp(x+dx).*sin(x+dx)) + ...
c_right(3).*(exp(x+2*dx).*sin(x+2*dx)) + c_right(4).*(exp(x+3*dx).*sin(x+3*dx)) + ...
c_right(5).*(exp(x+4*dx).*sin(x+4*dx)) + c_right(6).*(exp(x+5*dx).*sin(x+5*dx)))./(dx.^2);
% Absolute error between second order derivative of analytical and skewed right-sided scheme
skewed_right_error = abs(skewed_right - analytical_derivative);
% Left Sided Skewed
a_left = [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];
b_left = [0;0;1;0;0;0];
c_left = inv(a_left) * b_left;
% Calculating second order numerical derivative using Central differencing
skewed_left = (c_left(1).*(exp(x-5*dx).*sin(x-5*dx)) + c_left(2).*(exp(x-4*dx).*sin(x-4*dx)) + ...
c_left(3).*(exp(x-3*dx).*sin(x-3*dx)) + c_left(4).*(exp(x-2*dx).*sin(x-2*dx)) + ...
c_left(5).*(exp(x-dx).*sin(x-dx)) + c_left(6).*(exp(x)*sin(x)))./(dx.^2);
% Absolute error between second order derivative of analytical and skewed left-sided scheme
skewed_left_error = abs(skewed_left - analytical_derivative);
% Plotting the results
loglog(dx,central_error,'color','b','linewidth',2);
hold on
loglog(dx,skewed_right_error,'color','r','linewidth',2);
hold on
loglog(dx,skewed_left_error,'color','c','linewidth',2);
grid on
legend('Central Differencing','Right-Sided Skewed','Left-Sided Skewed')
xlabel('Range of dx')
ylabel('Errors')
title('Error Evaluation')
RESULT
CONCLUSION
- It has been concluded that the errors in central differencing method is much less than the other two schemes.
- The left-sided skewed and right sided skewed schemes shows nearly same amount of errors.
The use of skewed scheme is to compute the grids at the corners also which is not possible in central differencing scheme