Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Derivation of 4th order approximation of 2nd order derivative
required number of nodes = order of approximation +order of derivative-1
required number of nodes= 4+2-1=5
1. Central Difference Scheme: (five nodes are i-2, i-1, i, i+1, i+2)
`(del^2f)/(delx^2)=af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`----------(1)
now using Taylor series,
`af(i-2)=a[f(i)-2 Deltax f^'(i)+(2Deltax)^2/(2!)f^('')(i)-(2Deltax)^3/(3!)f^(''')(i)+(2Deltax)^4/(4!)f^('''')(i)-(2Deltax)^5/(5!)f^(''''')(i)+(2Deltax)^6/(6!)f^('''''')(i)]`
`bf(i-1)=b[f(i)- Deltax f^'(i)+(Deltax)^2/(2!)f^('')(i)-(Deltax)^3/(3!)f^(''')(i)+(Deltax)^4/(4!)f^('''')(i)-(Deltax)^5/(5!)f^(''''')(i)+(Deltax)^6/(6!)f^('''''')(i)]`
`cf(i)=cf(i)`
`df(i+1)=d[f(i)+ Deltax f^'(i)+(Deltax)^2/(2!)f^('')(i)+(Deltax)^3/(3!)f^(''')(i)+(Deltax)^4/(4!)f^('''')(i)+(Deltax)^5/(5!)f^(''''')(i)+(Deltax)^6/(6!)f^('''''')(i)]`
`ef(i+2)=e[f(i)+2 Deltax f^'(i)+(2Deltax)^2/(2!)f^('')(i)+(2Deltax)^3/(3!)f^(''')(i)+(2Deltax)^4/(4!)f^('''')(i)+(2Deltax)^5/(5!)f^(''''')(i)+(2Deltax)^6/(6!)f^('''''')(i)]`
| Table 1: Taylor table of "Central difference scheme" | |||||||
|
`f(i)` |
`Deltaxf^'(i)` |
`(Deltax)^2f^('')(i)` ` | `(Deltax)^3f^(''')(i)` | `(Deltax)^4f^('''')(i)` | `(Deltax)^5f^(''''')(i)` | `(Deltax)^6f^('''''')(i)` | |
| `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` | ||
| coefficients sum | 0 | 0 | 1 | 0 | 0 | ||
we are finding the second-order derivative`(del^2f)/(delx^2)`. so the coefficients of `f(i)`, `f^'(i)`,`f^(''')(i)`and `f^('''')(i)`should be equal to zero and coefficient of `f^('')(i)` must be equal to 1.
`a+b+c+d+e=0`
`-2a-b+0+d+2e=0`
`(4a)/2+b/2+0+d/2+(4e)/2=1`
`-(8a)/6-b/6+0+d/6+(8e)/6=0`
`(16a)/24+b/24+0+d/24+(16e)/24=0`
now convert these five equations into matrix form AX=B
`[[1 ,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]][[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`
now solve this matrix in MATLAB to calculate coefficients.
here a=-0.0833, b=1.3333, c=-2.5000, d=1.3333, f=-0.0833
4th order approximation of 2nd derivative by using central difference scheme
`(del^2f)/(delx^2)=(-0.0833f(i-2)+1.3333f(i-1)-2.5000f(i)+1.3333f(i+1)-0.0833f(i+2))/(Deltax^2)`
2. Skewed right-side difference scheme:
required number of nodes = order of approximation +order of derivative
required number of nodes= 4+2=6 (i.e i, i+1, i+2, i+3, i+4, i+5)
`(del^2f)/(delx^2)=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)`
from taylor series
`af(i)=af(i)`
`bf(i+1)=b[f(i)+ Deltax f^'(i)+(Deltax)^2/(2!)f^('')(i)+(Deltax)^3/(3!)f^(''')(i)+(Deltax)^4/(4!)f^('''')(i)+(Deltax)^5/(5!)f^(''''')(i)+(Deltax)^6/(6!)f^('''''')(i)]`
`cf(i+2)=c[f(i)+2 Deltax f^'(i)+(2Deltax)^2/(2!)f^('')(i)+(2Deltax)^3/(3!)f^(''')(i)+(2Deltax)^4/(4!)f^('''')(i)+(2Deltax)^5/(5!)f^(''''')(i)+(2Deltax)^6/(6!)f^('''''')(i)]`
`df(i+3)=d[f(i)+3 Deltax f^'(i)+(3Deltax)^2/(2!)f^('')(i)+(3Deltax)^3/(3!)f^(''')(i)+(3Deltax)^4/(4!)f^('''')(i)+(3Deltax)^5/(5!)f^(''''')(i)+(3Deltax)^6/(6!)f^('''''')(i)]`
`ef(i+4)=e[f(i)+4 Deltax f^'(i)+(4Deltax)^2/(2!)f^('')(i)+(4Deltax)^3/(3!)f^(''')(i)+(4Deltax)^4/(4!)f^('''')(i)+(4Deltax)^5/(5!)f^(''''')(i)+(4Deltax)^6/(6!)f^('''''')(i)]`
`gf(i+5)=g[f(i)+5 Deltax f^'(i)+(5Deltax)^2/(2!)f^('')(i)+(5Deltax)^3/(3!)f^(''')(i)+(5Deltax)^4/(4!)f^('''')(i)+(5Deltax)^5/(5!)f^(''''')(i)+(5Deltax)^6/(6!)f^('''''')(i)]`
| Table 2: Taylor table of "Skewed right-side difference scheme" | |||||||
|
`f(i)` |
`Deltaxf^'(i)` |
`(Deltax)^2f^('')(i)` ` | `(Deltax)^3f^(''')(i)` | `(Deltax)^4f^('''')(i)` | `(Deltax)^5f^(''''')(i)` | `(Deltax)^6f^('''''')(i)` | |
| `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` | |
| coefficients sum | 0 | 0 | 1 | 0 | 0 | ||
we are finding the second-order derivative`(del^2f)/(delx^2)`. so the coefficients of `f(i)`, `f^'(i)`,`f^(''')(i)`and `f^('''')(i)`should be equal to zero and coefficient of `f^('')(i)` must be equal to 1.
`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+(625e)/24=0`
`b/120+(32c)/120+(243d)/120+(1024e)/120+(3125g)/120=0`
now convert these six equations into matrix form AX=B
`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,2,9/2,8,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]]`
now solve this matrix in MATLAB to calculate coefficients.
here a= 3.7500, b=-12.8333, c=17.8333, d=-13.0000, e=5.0833, g=-0.8333
4th order approximation of 2nd derivative by skewed right side difference scheme
`(del^2f)/(delx^2)=(3.7500f(i)-12.8333f(i+1)+17.8333f(i+2)-13.0000f(i+3)+5.0833f(i+4)-0.8333f(i+5))/(Deltax^2)`
3. Skewed left-side difference scheme:
required number of nodes = order of approximation +order of derivative
required number of nodes= 4+2=6 (i.e i-5, i-4, i-3, i-2, i-1, i)
`(del^2f)/(delx^2)=af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)`
from taylor series
`af(i-5)=a[f(i)-5 Deltax f^'(i)+(5Deltax)^2/(2!)f^('')(i)-(5Deltax)^3/(3!)f^(''')(i)+(5Deltax)^4/(4!)f^('''')(i)-(5Deltax)^5/(5!)f^(''''')(i)+(5Deltax)^6/(6!)f^('''''')(i)]`
`bf(i-4)=b[f(i)-4 Deltax f^'(i)+(4Deltax)^2/(2!)f^('')(i)-(4Deltax)^3/(3!)f^(''')(i)+(4Deltax)^4/(4!)f^('''')(i)-(4Deltax)^5/(5!)f^(''''')(i)+(4Deltax)^6/(6!)f^('''''')(i)]`
`cf(i-3)=c[f(i)-3 Deltax f^'(i)+(3Deltax)^2/(2!)f^('')(i)-(3Deltax)^3/(3!)f^(''')(i)+(3Deltax)^4/(4!)f^('''')(i)-(3Deltax)^5/(5!)f^(''''')(i)+(3Deltax)^6/(6!)f^('''''')(i)]`
`df(i-2)=d[f(i)-2 Deltax f^'(i)+(2Deltax)^2/(2!)f^('')(i)-(2Deltax)^3/(3!)f^(''')(i)+(2Deltax)^4/(4!)f^('''')(i)-(2Deltax)^5/(5!)f^(''''')(i)+(2Deltax)^6/(6!)f^('''''')(i)]`
`ef(i-1)=e[f(i)- Deltax f^'(i)+(Deltax)^2/(2!)f^('')(i)-(Deltax)^3/(3!)f^(''')(i)+(Deltax)^4/(4!)f^('''')(i)-(Deltax)^5/(5!)f^(''''')(i)+(Deltax)^6/(6!)f^('''''')(i)]`
`gf(i)=gf(i)`
| Table 3: Taylor table of "Skewed left-side difference scheme" | |||||||
|
`f(i)` |
`Deltaxf^'(i)` |
`(Deltax)^2f^('')(i)` ` | `(Deltax)^3f^(''')(i)` | `(Deltax)^4f^('''')(i)` | `(Deltax)^5f^(''''')(i)` | `(Deltax)^6f^('''''')(i)` | |
| `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` | |
we are finding the second-order derivative`(del^2f)/(delx^2)`. so the coefficients of `f(i)`, `f^'(i)`,`f^(''')(i)`and `f^('''')(i)`should be equal to zero and coefficient of `f^('')(i)` must be equal to 1.
`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`
now convert these six equations into matrix form AX=B.
`[[1,1,1,1,1,1],[-5,-4,-3,-2,-1,0],[25/2,8,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],[g]]=[[0],[0],[1],[0],[0],[0]]`
now solve this matrix in MATLAB to calculate coefficients.
4th order approximation of 2nd order derivative by using skewed left side difference scheme
`(del^2f)/(delx^2)=(-0.8333f(i-5)+5.0833f(i-4)-13.0000f(i-3)+17.8333f(i-2)-12.8333f(i-1)+3.7500f(i))/(Deltax^2)`
clear all
close all
clc
% given function is f(x)=exp(x)*cos(x)
x=pi/3;
dx=linspace(pi/40,pi/4000,20);
for i=1:length(dx)
central_differencing(i)=central_difference(x,dx(i));
right_sided_skew(i)=skewed_right_side_difference(x,dx(i));
left_sided_skew(i)=skewed_left_side_difference(x,dx(i));
end
%% ploting
loglog(dx,central_differencing,'r')
hold on
plot(dx,right_sided_skew,'b')
hold on
plot(dx,left_sided_skew,'g')
xlabel('dx')
ylabel('Truncation error')
legend('central difference','skewed right side','skewed left side')
grid on
%% function of central difference scheme
function error_central_differencing=central_difference(x,dx)
analytic_derivative=-2*exp(x)*sin(x);
a=-0.8333;
b=1.3333;
c=2.50000;
d=1.3333;
e=0.8333;
central_differencing_fourth_order_approximation=((a*exp(x-2*dx)*cos(x-2*dx))+(b*exp(x-dx)*cos(x-dx))+(c*exp(x)*cos(x))+(d*exp(x+dx)*cos(x+dx))+(e*exp(x+2*dx)*cos(x+2*dx)))/dx^2;
error_central_differencing=abs(central_differencing_fourth_order_approximation-analytic_derivative);
end
%% function of skewed right side difference scheme
function error_skewed_right_side_difference=skewed_right_side_difference(x,dx)
analytic_derivative=-2*exp(x)*sin(x);
a=3.7500;
b=-12.8333;
c=17.8333;
d=-13.0000;
e=5.0833;
g=-0.8333;
skewed_right_side_fourth_order_approximation=((a*exp(x)*cos(x))+(b*exp(x+dx)*cos(x+dx))+(c*exp(x+2*dx)*cos(x+2*dx))+(d*exp(x+3*dx)*cos(x+3*dx))+(e*exp(x+4*dx)*cos(x+4*dx))+(g*exp(x+5*dx)*cos(x+5*dx)))/dx^2;
error_skewed_right_side_difference=abs(skewed_right_side_fourth_order_approximation-analytic_derivative);
end
%% function of skewed left side difference scheme
function error_skewed_left_side_difference=skewed_left_side_difference(x,dx)
analytic_derivative=-2*exp(x)*sin(x);
a=-0.8333;
b=-5.0833;
c=-13.0000;
d=-17.8333;
e=-12.8333;
g=3.7500;
skewed_left_side_fourth_order_approximation=((a*exp(x-5*dx)*cos(x-5*dx))+(b*exp(x-4*dx)*cos(x-4*dx))+(c*exp(x-3*dx)*cos(x-3*dx))+(d*exp(x-2*dx)*cos(x-2*dx))+(e*exp(x-dx)*cos(x-dx))+(g*exp(x)*cos(x)))/dx^2;
error_skewed_left_side_difference=abs(skewed_left_side_fourth_order_approximation-analytic_derivative);
end