Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

AIM: To write and execute a MATLAB program for 4th order approximations of the second order derivative using Taylor's Table method of the function exp(x)*cos(x) and compare the results with central difference approcimation and skewed schemes.

Deriving the 4th order approximations of the second-order derivative

1) Skewed right-sided difference

The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation

N= P+Q

where,

`N` is the number of nodes

`Q` is the order of derivative

`P` is the order of approximation

So, here we need 6 nodes in order to find the 2nd order derivative with the 4th order of approximation.

Consider the value of (del^2u)/(delx^2) with the 4th order of approximation is given by the equation. (del^2u)/(delx^2)=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5).

Now, using taylor series to find the values of discrete poits with repect to point `i`

`af(i)=af(i)`

`bf(i+1)=bf(i)+bDeltaxf'(i)+(bDeltax^2)/(2!)f''(i)+(bDeltax^3)/(3!)f'''(i)+(bDeltax^4)/(4!)f''''(i)+(bDeltax^5)/(5!)f'''''(i)+(bDeltax^6)/(6!)f''''''(i)+...`

`cf(i+2)=cf(i)+c(2Deltax)f'(i)+(c(2Deltax)^2)/(2!)f''(i)+(c(2Deltax)^3)/(3!)f'''(i)+(c(2Deltax)^4)/(4!)f''''(i)+(c(2Deltax)^5)/(5!)f'''''(i)+(c(2Deltax)^6)/(6!)f''''''(i)+...`

`df(i+3)=df(i)+d(3Deltax)f'(i)+(d(3Deltax)^2)/(2!)f''(i)+(d(3Deltax)^3)/(3!)f'''(i)+(d(3Deltax)^4)/(4!)f''''(i)+(d(3Deltax)^5)/(5!)f'''''(i)+(d(3Deltax)^6)/(6!)f''''''(i)+...`

`ef(i+4)=ef(i)+e(4Deltax)f'(i)+(e(4Deltax)^2)/(2!)f''(i)+(e(4Deltax)^3)/(3!)f'''(i)+(e(4Deltax)^4)/(4!)f''''(i)+(e(4Deltax)^5)/(5!)f'''''(i)+(e(4Deltax)^6)/(6!)f''''''(i)+...`

`gf(i+5)=gf(i)+g(5Deltax)f'(i)+g(5Deltax)^2/(2!)f''(i)+g(5Deltax)^3/(3!)f'''(i)+g(5Deltax)^4/(4!)f''''(i)+g(5Deltax)^5/(5!)f'''''(i)+g(5Deltax)^6/(6!)f''''''(i)+...`

Putting the LHS and RHS of the above equations we get

`af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)=f(i)[a+b+c+d+e+g]+Deltaxf'(i)[b+2c+3d+4e+5g]+Deltax^2f''(i)[b/2+(4c)/2+(9d)/2+(16e)/2+(25g)/2]+Deltax^3f'''(i)[b/6+(8c)/6+(27d)/6+(64e)/6+(125g)/6]+Deltax^4f''''(i)[b/24+(16c)/24+(81d)/24+(256e)/24+(625g)/24]+Deltax^5f'''''(i)[b/120+(32c)/120+(243d)/120+(1024e)/120+(3125g)/120]`

In the above equations, we want only the coefficients of `Deltax^2f''(i)` and we assume the coefficient to be one and for the remaining terms we assume the coefficients to be zero.

`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+(625g)/24=0`

`b/120+(32c)/120+(243d)/120+(1024e)/120+(3125g)/120=0`

TAYLOR TABLE

Function	`f(i)`	`Deltaf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f'''(i)`	`Deltax^5'''''(i)`	`Deltax^6''''''(i)`
`af(i)`	`a`	`0`	`0`	`0`	`0`	`0`	`0`
`bf(i+1)`	`b`	`b`	`b/2`	`b/6`	`b/24`	`b/120`	`b/720`
`cf(i+2)`	`c`	`2c`	`(4c)/2`	`(8c)/6`	`(16c)/24`	`(32c)/120`	`(64c)/720`
`df(i+3)`	`d`	`3d`	`(9d/2)`	`(27d)/6`	`(81d)/24`	`(243d)/120`	`(729d)/720`
`ef(i+4)`	`e`	`4e`	`(16e)/2`	`(64e)/6`	`(256e)/24`	`(1024e)/120`	`(4096e)/720`
`gf(i+5)`	`g`	`5g`	`(25g)/2`	`(125g)/6`	`(625g)/24`	`(3125g)/120`	`(15625g)/720`
Sum of Coefficient	0	0	1	0	0	0	?

the linear equation we get from the above table are solved using matrix inversion method,

Writing the coefficients in the form of a matrix we get,

`[[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]]`

Solving using MATLAB to find the values of the coefficients and we get,

2) Skewed left-sided difference

The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation

`N=P+Q`

where

`N` is the number of nodes

`Q` is the order of derivative

`P` is the order of approximation

So, here we need 6 nodes in order to find the 2nd order derivative with the 4th order of approximation.

Consider the value of `(del^2y)/(delx^2)` with the 4th order of approximation is given by the equation

`(del^2y)/(delx^2)=af(i)+bf(i-1)+cf(i-2)+df(i-3)+ef(i-4)+gf(i-5)`

TAYLOR TABLE

Function	`f(i)`	`Deltaf'(i)`	`Deltax^2f''(i)`	`Deltax3f'''(i)`	`Deltax4f''''(i)`	`Deltax^5'''''(i)`	`Deltax^6''''''(i)`
`af(i)`	`a`	`0`	`0`	`0`	`0`	`0`	`0`
`bf(i-1)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`	`b/720`
`cf(i-2)`	`c`	`-2c`	`(4c)/2`	`-(8c)/6`	`(16c)/24`	`-(32c)/120`	`(64c)/720`
`df(i-3)`	`d`	`-3d`	`(9d)/2`	`-(27d)/6`	`(81d)/24`	`-(243d)/120`	`(729d)/720`
`ef(i-4)`	`e`	`-4e`	`(16e)/2`	`-(64e)/6`	`(256e)/24`	`-(1024e)/120`	`(4096e)720`
`gf(i-5)`	`g`	`-5g`	`(25g)/2`	`-(125g)/6`	`(625g)/24`	`-(3125g)/120`	`(15625g)/720`
Sum of Coefficient	0	0	1	0	0	0	?

the linear equation we get from the above table are solved using matrix inversion method,

Writing the coefficients in the form of a matrix we get,

`[[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]]`

Solving using MATLAB to find the values of the coefficients and we get,

3. Central Difference

The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation

`N=P+Q-1`

where

`N` is the number of nodes

`Q` is the order of derivative

`P` is the order of approximation

So, here we need 5 nodes in order to find the 2nd order derivative with the 4th order of approximation.

Consider the value of `(del^2y)/(delx^2)` with the 4th order of approximation is given by the equation

`(del^2y)/(delx^2)=af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`

Function	`f(i)`	`Deltaf'(i)`	`Deltax^2f''(i)`	`Deltax3f'''(i)`	`Deltax4f''''(i)`	`Deltax^5'''''(i)`
`af(i-2)`	`a`	`-2a`	`(4a)/2`	`-(8a)/6`	`(16a)/24`	`-(32a)/120`
`bf(i-1)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`
`cf(i)`	`c`	`0`	`0`	`0`	`0`	`0`
`df(i+1)`	`d`	`d`	`d/2`	`d/6`	`d/24`	`-c/120`
`ef(i+2)`	`e`	`2e`	`(4e)/2`	`-(8e)/6`	`(16e)/24`	`(32e)/120`
Sum of Coefficient	0	0	1	0	0	?

the linear equation we get from the above table are solved using matrix inversion method,

Writing the coefficients in the form of a matrix we get,

`[[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],[-32/120, -1/120, 0, 1/120, 32/120]][[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

Solving using MATLAB to find the values of the coefficients and we get,

MATLAB CODE:

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,'r')
hold on
loglog(dx,right_sided_skew,'g')
hold on
loglog(dx,left_sided_skew,'b')
xlabel('range of dx')
ylabel('truncation error')
title('Fourth order approximation')
legend('central differencing','skewed right sided difference','skewed left sided difference')

CODE FOR CENTRAL DIFFERENCING:

function Error_central_differencing = fourth_order_central_differencing(x,dx)

% analytic function = exp(x)*cos(x)

% analytic 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
Error_central_differencing = abs(central_differencing- analytical_derivative);
end

CODE FOR SKEWED RIGHT SIDED DIFFERENCE:

function Error_right_sided_skew = fourth_order_right_sided_skew(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

analytical_derivative = -2*exp(x)*sin(x);

%numerical derivative

%central differencing
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
Error_right_sided_skew = abs(right_sided_skew- analytical_derivative);
end

CODE FOR SKEWED LEFT SIDED DIFFERENCE:

function Error_left_sided_skew = fourth_order_left_sided_skew(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

analytical_derivative = -2*exp(x)*sin(x);

%numerical derivative

%central differencing
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
Error_left_sided_skew = abs(left_sided_skew- analytical_derivative);
end

Output:

Why Skewed schemes are useful?

When we dont have any nodes at the boundary then skewed schemes are useful to approximate to improve the accuracy but when we have internal nodes available then Centered ones are good though.

Conclusion:

The comparison of absolute error for the left, right and central differencing scheme can be seen in the below graph. The graph shows that the central scheme is most accurate compared to the skewed differencing. The curves also show that when the value of dx is very small or large the approximating function becomes unstable. This instability is caused due to the combination of various errors such as roundoff error, residual error and truncation error.

From the below graph it is seen that the range of dx where the function is stable ranges from 0.1 to 0.01, beyond this range we see an abrupt rise or fall.