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 using the Taylor table method.
- Central difference.
- Skewed right-sided difference.
- Skewed left-sided difference.
OBJECTIVES:
- Prove that obtained skewed schemes are fourth-order accurate.
- 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.
- The plot compares the absolute error between the above-mentioned schemes.
DERIVATION OF 4th ORDER APPROXIMATION OF A 2nd ORDER DERIVATIVE USING TAYLOR TABLE METHOD
Using Skewed right-sided difference
In order to derive the 4th order skewness for the right-hand side from point (i), we have to consider five points on the right of (i). Let us consider (i+1), (i+2), (i+3), (i+4), and (i+5) with spacing of `Delta x`.
We need 4th order approximation hence we take stencils, hence the hypothesis equation we get
`(del^2u)/(delx^2)=a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)`
From the above hypothesis equation to find (a,b,c,d,e,g) we are using the Taylor table method
Firstly table we do the Taylor series expansion for f(i+1), f(i+2), f(i+3), f(i+4), and f(i+5)
`a*f(i)=a*f(i)`
`b*f(i+1)= b*f(i)+ ((b)*f^1(i)*(Delta x))/(1!)+((b)*f^2(i)*(Delta x)^2)/(2!)+((b)*f^3(i)*(Delta x)^3)/(3!)+((b)*f^4(i)*(Delta x)^4)/(4!)`
`+((b)*f^5(i)*(Delta x)^5)/(5!)+((b)*f^6(i)*(Delta x)^6)/(6!)`
`c*f(i+2)= c*f(i)+ ((c)*f^1(i)*(2*Delta x))/(1!)+((c)*f^2(i)*(2*Delta x)^2)/(2!)+((c)*f^3(i)*(2*Delta x)^3)/(3!)+((c)*f^4(i)*(2*Delta x)^4)/(4!)`
`+((c)*f^5(i)*(2*Delta x)^5)/(5!)+((c)*f^6(i)*(2*Delta x)^6)/(6!)`
`d*f(i+3)= d*f(i)+ ((d)*f^1(i)*(3*Delta x))/(1!)+((d)*f^2(i)*(3*Delta x)^2)/(2!)+((d)*f^3(i)*(3*Delta x)^3)/(3!)+((d)*f^4(i)*(3*Delta x)^4)/(4!)`
`+((d)*f^5(i)*(3*Delta x)^5)/(5!)+((d)*f^6(i)*(3*Delta x)^6)/(6!)`
`e*f(i+4)= e*f(i)+ ((e)*f^1(i)*(4*Delta x))/(1!)+((e)*f^2(i)*(4*Delta x)^2)/(2!)+((e)*f^3(i)*(4*Delta x)^3)/(3!)+((e)*f^4(i)*(4*Delta x)^4)/(4!)`
`+((e)*f^5(i)*(4*Delta x)^5)/(5!)+((e)*f^6(i)*(4*Delta x)^6)/(6!)`
`g*f(i+5)= g*f(i)+ ((g)*f^1(i)*(5*Delta x))/(1!)+((g)*f^2(i)*(5*Delta x)^2)/(2!)+((g)*f^3(i)*(5*Delta x)^3)/(3!)+((g)*f^4(i)*(5*Delta x)^4)/(4!)`
`+((g)*f^5(i)*(5*Delta x)^5)/(5!)+((g)*f^6(i)*(5*Delta x)^6)/(6!)`
By arranging the Taylor series coefficients in the Taylor table as shown below.
| `f(i)` | `Deltax*f(i)` | `Deltax^2*f^2(i)` | `Deltax^3*f^3(i)` | `Deltax^4*f^4(i)` | `Deltax^5*f^5(i)` | `Deltax^6*f^6(i)` | |
| `a*f(i)` | a | 0 | 0 | 0 | 0 | 0 | 0 |
| `b*f(i+1)` | b | b | `b/2` | `b/6` | `b/24` | `b/120` | `b/720` |
| `c*f(i+2)` | c | `2*c` | `(4*c)/2` | `(8*c)/6` | `(16*c)/24` | `(32*c)/120` | `(64*c)/720` |
| `d*f(i+3)` | d | `3*d` | `(9*d)/2` | `(27*d)/6` | `(81*d)/24` | `(243*d)/120` | `(729*d)/720` |
| `e*f(i+4)` | e | `4*e` | `(16*e)/2` | `(64*e)/6` | `(256*e)/24` | `(1024*e)/120` | `(4096*e)/720` |
| `g*f(i+5)` | g | `5*g` | `(25*g)/2` | `(125*g)/6` | `(625*g)/24` | `(3125*g)/120` | `(15625*g)/720` |
| 0 | 0 | 1 | 0 | 0 | 0 |
By adding the above coefficients from the tables we get the following equations as shown below.
`a+b+c+d+e+g=0`
`0+b+2*c+3*d+4*e+5*g=0`
`0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2=1`
`0+b/6+(8*c)/6+(27*d)/6+(64*e)/6+(125*g)/6=0`
`0+b/24+(16*c)/24+(81*d)/24+(256*e)/24+(625*g)/24=0`
`0+b/120+(32*c)/120+(243*d)/120+(1024*e)/120+(3125*g)/120=0`
Here we have six unknown (a,b,c,d,e,g) coefficients, to get the solutions for the above equations we have to solve the above in matrix form AX=B
`[[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]]`
After solving the above matrix using the MATLAB we get all coefficients as shown below
MATLAB CODE
A=[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=[0 0 1 0 0 0];
X=B.*inv(A)`[[a=3.75],[b=-12.8333],[c=17.8333],[d=-13.0],[e=5.0833],[g=-0.8333]]`
From Taylor's table(LHS),
`(del^2u)/(delx^2)=a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5)`for this hypothesis equation we are dividing `Delta x^2` we get
`(del^2u)/(delx^2)=(a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+g*f(i+5))/(Delta x^2)` this RHS of the skewness right hand difference. By substituting all coefficient values we get
`(del^2u)/(delx^2)=(3.75*f(i)-12.8333*f(i+1)+17.8333c*f(i+2)-13*f(i+3)+5.0833*f(i+4)-0.8333*f(i+5))/(Delta x^2)`
From Taylor's table(RHS),
`(del^2u)/(delx^2)=(a+b+c+d+e+g)*f(i)+(0+b+2*c+3*d+4*e+5*g)*Delta x*f(i)+(0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2)*Delta x^2*f^2(i)`
`+(0+b/6+(8*c)/6+(27*d)/6+(64*e)/6+(125*g)/6)*Delta x^3*f^3(i)+(0+b/24+(16*c)/24+(81*d)/24+(256*e)/24+(625*g)/24)*Delta x^4*f^4(i)`
`+(0+b/120+(32*c)/120+(243*d)/120+(1024*e)/120+(3125*g)/120)*Deltax^5*f^5(i)`
To get the second order derivative for the RHS side of skewness right-hand difference for `f^2(i)`is equal to 1, for `f(i)`, `f^1(i)`, `f^3(i)`, `f^4(i)`, `f^5(i)` equal to 0 by substituting the condition to the above equation we get
`(del^2u)/(delx^2)=(0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2)*Delta x^2*f^2(i)+(0+b/720+(64*c)/720+(729*d)/720+(4096*e)/720+(15625*g)/720)*Deltax^6*f^6(i)`
As discussed before by dividing `Delta x^2` we get
`(del^2u)/(delx^2)=(0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2)*f^2(i)+(0+b/720+(64*c)/720+(729*d)/720+(4096*e)/720+(15625*g)/720)*Deltax^4*f^4(i)`
Hence `Delta x^4` denotes the truncation error. The order of truncation error is 4. we know that the order of approximation is defined by the truncation error term. Hence we can say that the skewness right-hand difference scheme is 4th order accurate.
Using Skewed left-sided difference
In order to derive the 4th order skewness for the left-hand side from point (i), we have to consider five points on the right of (i). Let us consider (i-1), (i-2), (i-3), (i-4), and (i-5) with spacing of `Delta x`.
We need 4th order approximation hence we take stencils, hence the hypothesis equation we get
`(del^2u)/(delx^2)=a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5)`
From the above hypothesis equation to find (a,b,c,d,e,g) we are using the Taylor table method
Firstly table we do the Taylor series expansion for f(i-1), f(i-2), f(i-3), f(i-4), and f(i-5)
`a*f(i)=a*f(i)`
`b*f(i-1)= b*f(i)- ((b)*f^1(i)*(Delta x))/(1!)+((b)*f^2(i)*(Delta x)^2)/(2!)-((b)*f^3(i)*(Delta x)^3)/(3!)+((b)*f^4(i)*(Delta x)^4)/(4!)`
`-((b)*f^5(i)*(Delta x)^5)/(5!)+((b)*f^6(i)*(Delta x)^6)/(6!)`
`c*f(i-2)= c*f(i)- ((c)*f^1(i)*(2*Delta x))/(1!)+((c)*f^2(i)*(2*Delta x)^2)/(2!)-((c)*f^3(i)*(2*Delta x)^3)/(3!)+((c)*f^4(i)*(2*Delta x)^4)/(4!)`
`-((c)*f^5(i)*(2*Delta x)^5)/(5!)+((c)*f^6(i)*(2*Delta x)^6)/(6!)`
`d*f(i-3)= d*f(i)- ((d)*f^1(i)*(3*Delta x))/(1!)+((d)*f^2(i)*(3*Delta x)^2)/(2!)-((d)*f^3(i)*(3*Delta x)^3)/(3!)+((d)*f^4(i)*(3*Delta x)^4)/(4!)`
`-((d)*f^5(i)*(3*Delta x)^5)/(5!)+((d)*f^6(i)*(3*Delta x)^6)/(6!)`
`e*f(i-4)= e*f(i)- ((e)*f^1(i)*(4*Delta x))/(1!)+((e)*f^2(i)*(4*Delta x)^2)/(2!)-((e)*f^3(i)*(4*Delta x)^3)/(3!)+((e)*f^4(i)*(4*Delta x)^4)/(4!)`
`-((e)*f^5(i)*(4*Delta x)^5)/(5!)+((e)*f^6(i)*(4*Delta x)^6)/(6!)`
`g*f(i-5)= g*f(i)- ((g)*f^1(i)*(5*Delta x))/(1!)+((g)*f^2(i)*(5*Delta x)^2)/(2!)-((g)*f^3(i)*(5*Delta x)^3)/(3!)+((g)*f^4(i)*(5*Delta x)^4)/(4!)`
`-((g)*f^5(i)*(5*Delta x)^5)/(5!)+((g)*f^6(i)*(5*Delta x)^6)/(6!)`
By arranging the Taylor series coefficients in the Taylor table as shown below.
| `f(i)` | `Deltax*f(i)` | `Deltax^2*f^2(i)` | `Deltax^3*f^3(i)` | `Deltax^4*f^4(i)` | `Deltax^5*f^5(i)` | `Deltax^6*f^6(i)` | |
| `a*f(i)` | a | 0 | 0 | 0 | 0 | 0 | 0 |
| `b*f(i-1)` | b | -b | `b/2` | `-b/6` | `b/24` | `-b/120` | `b/720` |
| `c*f(i-2)` | c | `-2*c` | `(4*c)/2` | `(-8*c)/6` | `(16*c)/24` | `(-32*c)/120` | `(64*c)/720` |
| `d*f(i-3)` | d | `-3*d` | `(9*d)/2` | `(-27*d)/6` | `(81*d)/24` | `(-243*d)/120` | `(729*d)/720` |
| `e*f(i-4)` | e | `-4*e` | `(16*e)/2` | `(-64*e)/6` | `(256*e)/24` | `(-1024*e)/120` | `(4096*e)/720` |
| `g*f(i-5)` | g | `-5*g` | `(25*g)/2` | `(-125*g)/6` | `(625*g)/24` | `(-3125*g)/120` | `(15625*g)/720` |
| 0 | 0 | 1 | 0 | 0 | 0 |
By adding the above coefficients from the tables we get the following equations as shown below.
`a+b+c+d+e+g=0`
`0-b-2*c-3*d-4*e-5*g=0`
`0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2=1`
`0-b/6-(8*c)/6-(27*d)/6-(64*e)/6-(125*g)/6=0`
`0+b/24+(16*c)/24+(81*d)/24+(256*e)/24+(625*g)/24=0`
`0-b/120-(32*c)/120-(243*d)/120-(1024*e)/120-(3125*g)/120=0`
Here we have six unknown (a,b,c,d,e,g) coefficients, to get the solutions for the above equations we have to solve the above in matrix form AX=B
`[[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]]`
After solving the above matrix using the MATLAB we get all coefficients as shown below
MATLAB CODE
A=[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=[0 0 1 0 0 0];
X=B.*inv(A)`[[a=3.75],[b=-12.8333],[c=17.8333],[d=-13.0],[e=5.0833],[g=-0.8333]]`
From Taylor's table(LHS),
`(del^2u)/(delx^2)=a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5)`for this hypothesis equation we are dividing `Delta x^2` we get
`(del^2u)/(delx^2)=(a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+g*f(i-5))/(Delta x^2)` this RHS of the skewness left hand difference. By substituting all coefficient values we get
`(del^2u)/(delx^2)=(3.75*f(i)-12.8333*f(i-1)+17.8333c*f(i-2)-13*f(i-3)+5.0833*f(i-4)-0.8333*f(i-5))/(Delta x^2)`
From Taylor's table(RHS),
`(del^2u)/(delx^2)=(a+b+c+d+e+g)*f(i)+(0-b-2*c-3*d-4*e-5*g)*Delta x*f(i)+(0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2)*Delta x^2*f^2(i)`
`+(0-b/6-(8*c)/6-(27*d)/6-(64*e)/6-(125*g)/6)*Delta x^3*f^3(i)+(0+b/24+(16*c)/24+(81*d)/24+(256*e)/24+(625*g)/24)*Delta x^4*f^4(i)`
`+(0-b/120-(32*c)/120-(243*d)/120-(1024*e)/120-(3125*g)/120)*Deltax^5*f^5(i)`
To get the second order derivative for the RHS side of skewness left-hand difference for `f^2(i)`is equal to 1, for `f(i)`, `f^1(i)`, `f^3(i)`, `f^4(i)`, `f^5(i)` equal to 0 by substituting the condition to the above equation we get
`(del^2u)/(delx^2)=(0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2)*Delta x^2*f^2(i)+(0+b/720+(64*c)/720+(729*d)/720+(4096*e)/720+(15625*g)/720)*Deltax^6*f^6(i)`
As discussed before by dividing `Delta x^2` we get
`(del^2u)/(delx^2)=(0+b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2)*f^2(i)+(0+b/720+(64*c)/720+(729*d)/720+(4096*e)/720+(15625*g)/720)*Deltax^4*f^4(i)`
Hence `Delta x^4` denotes the truncation error. The order of truncation error is 4. we know that the order of approximation is defined by the truncation error term. Hence we can say that the skewness left-hand difference scheme is 4th order accurate.
Using a central difference scheme
For computing the 4th order approximation at point(i), we have to take the information from the four points surrounding that we are going to take two points on the left of (i) and two points on the right of (i). let the points be (i-2), (i-1), (i+1), (i+2) and there distances are`Delta x`.
We need 4th order approximation hence we take stencils, hence the hypothesis equation we get
`(del^2u)/(delx^2)=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)`
From the above hypothesis equation to find (a,b,c,d,e,g) we are using the Taylor table method
Firstly table we do the Taylor series expansion for f(i-2), f(i-1), f(i), f(i+1), and f(i+2)
`a*f(i-2)= a*f(i)- ((a)*f^1(i)*(2*Delta x))/(1!)+((a)*f^2(i)*(2*Delta x)^2)/(2!)-((a)*f^3(i)*(2*Delta x)^3)/(3!)+((a)*f^4(i)*(2*Delta x)^4)/(4!)`
`-((a)*f^5(i)*(2*Delta x)^5)/(5!)+((a)*f^6(i)*(2*Delta x)^6)/(6!)`
`b*f(i-1)= b*f(i)- ((b)*f^1(i)*(Delta x))/(1!)+((b)*f^2(i)*(Delta x)^2)/(2!)-((b)*f^3(i)*(Delta x)^3)/(3!)+((b)*f^4(i)*(Delta x)^4)/(4!)`
`-((b)*f^5(i)*(Delta x)^5)/(5!)+((b)*f^6(i)*(Delta x)^6)/(6!)`
`c*f(i)=c*f(i)`
`d*f(i+1)= d*f(i)+ ((d)*f^1(i)*(Delta x))/(1!)+((d)*f^2(i)*(Delta x)^2)/(2!)+((d)*f^3(i)*(Delta x)^3)/(3!)+((d)*f^4(i)*(Delta x)^4)/(4!)`
`+((d)*f^5(i)*(Delta x)^5)/(5!)+((d)*f^6(i)*(Delta x)^6)/(6!)`
`e*f(i+2)= e*f(i)+ ((e)*f^1(i)*(2*Delta x))/(1!)+((e)*f^2(i)*(2*Delta x)^2)/(2!)+((e)*f^3(i)*(2*Delta x)^3)/(3!)+((e)*f^4(i)*(2*Delta x)^4)/(4!)`
`+((e)*f^5(i)*(2*Delta x)^5)/(5!)+((e)*f^6(i)*(2*Delta x)^6)/(6!)`
By arranging the Taylor series coefficients in the Taylor table as shown below.
| `f(i)` | `Deltax*f(i)` | `Deltax^2*f^2(i)` | `Deltax^3*f^3(i)` | `Deltax^4*f^4(i)` | `Deltax^5*f^5(i)` | `Deltax^6*f^6(i)` | |
| `a*f(i-2)` | a | `-2*a` | `(4*a)/2` | `-(8*a)/6` | `(16*a)/24` | `-(32*a)/120` | `(64*a)/240` |
| `b*f(i-1)` | b | -b | `b/2` | `-b/6` | `b/24` | `-b/120` | `b/240` |
| `c*f(i)` | c | 0 | 0 | 0 | 0 | 0 | 0 |
| `d*f(i+1)` | d | d | `d/2` | `d/6` | `d/24` | `d/120` | `d/240` |
| `e*f(i+2)` | e | `2*e` | `(4*e)/2` | `(8*e)/6` | `(16*e)/24` | `(32*e)/120` | `(64*e)/240` |
| 0 | 0 | 1 | 0 | 0 | 0 |
By adding the above coefficients from the tables we get the following equations as shown below.
`a+b+c+d+e+g=0`
`-2*a-b+0+d+2*e=0`
`2*a+b/2+0+d/2+e*2=1`
`-(8*a)/6-b/6+0+d/6+(8*e)/6=0`
`(16*a)/24+b/24+0+d/24+(16*e)/24=0`
Here we have six unknown (a,b,c,d,e) coefficients, to get the solutions for the above equations we have to solve the above in 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]]`
After solving the above matrix using the MATLAB we get all coefficients as shown below
MATLAB CODE
A=[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];
B=[0 0 1 0 0 ];
X=B.*inv(A)`[[a=-0.0833],[b=1.333],[c=-2.5],[d=1.3333],[e=-0.0833]]`
From Taylor's table(LHS),
`(del^2u)/(delx^2)=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)`for this hypothesis equation we are dividing `Delta x^2` we get
`(del^2u)/(delx^2)=(a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2))/(Delta x^2)` this RHS of the central difference. By substituting all coefficient values we get
`(del^2u)/(delx^2)=(-0.0833*f(i-2)+1.3333*f(i-1)-2.5*f(i)+1.3333*f(i+1)-0.0833*f(i+2))/(Delta x^2)`
From Taylor's table(RHS),
`(del^2u)/(delx^2)=(a+b+c+d+e)*f(i)+(-2*a-b+0+d+2*e)*Delta x*f(i)+(2*a+b/2+0+d/2+e*2)*Delta x^2*f^2(i)`
`+(-(8*a)/6-b/6+0+d/6+(8*e)/6)*Delta x^3*f^3(i)+((16*a)/24+b/24+0+d/24+(16*e)/24)*Delta x^4*f^4(i)+((-32*a)/120-b/120+0+d/120+(32*e)/120)*Delta x^5*f^5)+(62*a)/240+b/240+0+d/240+(62*e/240)*Delta x^6*f^6)`
To get the second order derivative for the RHS side of central difference for `f^2(i)`is equal to 1, for `f(i)`, `f^1(i)`, `f^3(i)`, `f^4(i)`, `f^5(i)` equal to 0 by substituting the condition to the above equation we get
`(del^2u)/(delx^2)=(2*a+b/2+0+d/2+e*2)*Delta x^2*f^2(i)+(62*a)/240+b/240+0+d/240+(62*e/240)*Delta x^6*f^6)`
As discussed before by dividing `Delta x^2` we get
`(del^2u)/(delx^2)=(2*a+b/2+0+d/2+e*2)*f^2(i)+(62*a)/240+b/240+0+d/240+(62*e/240)*Delta x^4*f^4)`
Hence `Delta x^4` denotes the truncation error. The order of truncation error is 4. we know that the order of approximation is defined by the truncation error term. Hence we can say that the central difference scheme is 4th order accurate.
MATLAB 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
RIGHT SKEWED DIFFERENCING
% Creating function for the 4th order approximation of the second-order derivative using skewed right-sided
function error_skewed_right_sided=skewed_right_sided(x,dx)
% analytical function f(x)=exp(x)*cos(x);
%f''(x)=-2*exp(x)*sin(x);
analytical_derivative=-2*exp(x)*sin(x);
%values of coefficients calculated from our Taylor table;
a=3.75000;
b=-12.83333;
c=17.83333;
d=-13.00000;
e=5.08333;
g=-0.83333;
%numerical derivative;
%skewed right-sided differencing=((a*f(x))+(b*f(x+dx))+(c*f(x+2*dx))+(d*f(x+3*dx))+(e*f(x+4*dx))+(g*f(x+5*dx)))/dx^2;
%skewed right-sided fourth order approximation;
skewed_right_sided=((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_sided=abs(skewed_right_sided-analytical_derivative);
endLEFT SKEWED DIFFERENCING
% Creating function for the 4th order approximation of the second-order derivative using skewed left-sided
function error_skewed_left_sided=skewed_left_sided(x,dx)
% analytical function f(x)=exp(x)*cos(x);
%f''(x)=-2*exp(x)*sin(x);
analytical_derivative=-2*exp(x)*sin(x);
%values of coefficients calculated from our Taylor table;
a=-0.83333;
b=5.08333;
c=-13.00000;
d=17.83333;
e=-12.83333;
g=3.75000;
%numerical derivative;
%skewed left-sided differencing=((a*f(x-5*dx))+(b*f(x-4*dx))+(c*f(x-3*dx))+(d*f(x-2*dx))+(e*f(x-dx))+(g*f(x)))/dx^2;
%skewed left-sided fourth order approximation;
skewed_left_sided=((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_sided=abs(skewed_left_sided-analytical_derivative);
endCENTRAL DIFFERENCING
% Creating function for the 4th order approximation of the second-order; derivative using central difference
function error_central_differencing = central_difference(x,dx)
% Analytical function f(x) = exp(x)*cos(x)
% f''(x) = -2*exp(x)*sin(x);
analytical_derivative = -2*exp(x)*sin(x);
% values of coefficients calculated from our Taylor table
a=-0.083333;
b=1.333333;
c=-2.500000;
d=1.333333;
e=-0.083333;
% numerical derivative
% central differencing =
% ((a*f(x-2*dx))+(b*f(x-dx))+(c*f(x))+(d*f(x+dx))+(e*f(x+2*dx)))/dx^2
%central difference fourth order approximation
central_differencing = ((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-analytical_derivative);
endMAIN PROGRAM
%calculating 4th order error of second-order derivative using central, skewed right-sided, and left-sided ;
clear all
close all
clc
x=pi/3;
dx=linspace(pi/4,pi/40000,40);
%creating 'for' loop;
for i=1:length(dx)
central_difference_error(i)=central_difference(x,dx(i));
skewed_right_side_error(i)=skewed_right_sided(x,dx(i));
skewed_left_side_error(i)=skewed_left_sided(x,dx(i));
end
%plotting graphs for dx vs errors for central, skewed right-sided, & skewed left-sided approximations
figure(1)
plot(dx,central_difference_error,'linewidth',2);
hold on
grid on
plot(dx,skewed_right_side_error,'r','linewidth',2);
hold on
plot(dx,skewed_left_side_error,'g','linewidth',2);
xlabel('dx');
ylabel('error');
legend('Central Difference Error','Skewed Right-Sided Error','Skewed Left-Sided Error');
figure(2)
loglog(dx,central_difference_error,'linewidth',2);
hold on
grid on
loglog(dx,skewed_right_side_error,'r','linewidth',2);
hold on
loglog(dx,skewed_left_side_error,'g','linewidth',2);
xlabel('dx');
ylabel('error');
legend('Central Difference Error','Skewed Right-Sided Error','Skewed Left-Sided Error');CODE EXPLANATION:
- Firstly we’re going to write programs to develop functions for central, skewed right-sided, and skewed left-sided errors at a particular point (x, dx), where dx represents spacing in between each node.
- These functions will help us to calculate errors over the range of ‘dx’. With the help of code ‘function out we’re going to make a function that we save in the directory and call according to our need in the main program (make sure to save the main program in the same directory in which your function is saved).
-
After creating the above functions, we save our main program in the same directory in which we have saved our three functions; central_difference, skewed_right_sided, and skewed_left_sided.
-
For now, we’re going to create an array for the range of ‘dx’. We’ll use ‘for loop’ over the range of ‘dx’ and use the functions central, skewed right-sided, and skewed left-sided difference 4th order error at a given point (x, dx).
-
we’ll get the errors for a range of ‘dx’ and also we’re going to plot the graph between ‘error’ vs. ‘dx’ for central, skewed right-sided, and skewed left-sided difference 4th order error.
why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?
As we can see in this problem central differencing plays a good role in reducing the error compared to the skewed scheme by taking the values of both sides but we cannot always apply the central differencing, especially at end node points or initial boundary conditions. In order to solve node points or initial boundary conditions left and right skewed plays a major role in determining the accurate solution for a given problem statement. Hence both scheme plays a vital role in determining the solution to a given problem statement.
RESULTS:
CONCLUSIONS:
Analyzing the graph between ‘error’ vs. ‘dx’ for central, skewed right-sided, and skewed left-sided difference 4th order error. We can easily notice that central difference approximation produces more accurate solutions. Hence produces lesser errors as compared to right-skewed and left-skewed difference schemes.