Taylor Table Method For Numerical Approximation Methods

Objective :

1. To find the fourth order approximations of second order derivative , using central and skewed scheme.
2. Compare all schemes.

Central Symmetric Scheme : Fourth Order Approximation of second order derivative

Since , Fourth order Approximation we need four discrete neighbour points .

Central , so same number of nodes on both sides.

Method:

`frac{d^2f}{dx^2}= a*f(i-2)+ b*f(i-1)+c*f(i) + d*f(i+1) + e*f(i+2)+...` ------->(1)

Using Taylor Series Expansion:

`a*f(i-2) approx a*f(x_i) -a*f^1(x_i)*2*h + a^2*frac{f^2(x_i)}{2!}*(2*h)^2 - a*frac{f^3(x_i)}{3!}*(2*h)^3 +a*frac{f^4(x_i)}{4!}*(2*h)^4 - a*frac{f^5(x_i)}{5!}*(2*h)^5+...` -->(2)

`b*f(i-1) approx b*f(x_i) -b*f^1(x_i)*h + b*frac{f^2(x_i)}{2!}*(*h)^2 - b*frac{f^3(x_i)}{3!}*(h)^3 +b*frac{f^4(x_i)}{4!}*(h)^4 - b*frac{f^5(x_i)}{5!}*(h)^5+...` -->(3)

`c*f(i) =c*f(x_i)` -->(4)

`b*f(i+1) approx b*f(x_i) +b*f^1(x_i)*h + b*frac{f^2(x_i)}{2!}*(*h)^2 + b*frac{f^3(x_i)}{3!}*(h)^3 +b*frac{f^4(x_i)}{4!}*(h)^4 + b*frac{f^5(x_i)}{5!}*(h)^5+...` -->(5)

`a*f(i+2) approx a*f(x_i) +a*f^1(x_i)*2*h + a^2*frac{f^2(x_i)}{2!}*(2*h)^2 + a*frac{f^3(x_i)}{3!}*(2*h)^3 +a*frac{f^4(x_i)}{4!}*(2*h)^4 + a*frac{f^5(x_i)}{5!}*(2*h)^5+...` -->(6)

Here taken , h = grid space.

Taking common terms and equating to L.H.S , represented in table

	`f(x_i)`	`f^1(x_i)*Deltax`	`f^2(x_i)*(Deltax)^2`	`f^3(x_i)(Deltax)^3`	`f^4(x_i)*(Deltax)^4`	`f^5(x_i)*(Deltax)^5`	`f^6(x_i)*(Deltax)^6`
`a*f_(i-2)`	a	`-2*a`	`2*a`	`frac{-8*a}{3}`	`frac{16*a}{24}`	`frac{-32*a}{120}`	`frac{2^6}{6!}*e`
`b*f_(i-1)`	b	`-b`	`frac{b}{2}`	`frac{-b}{6}`	`frac{b}{24}`	`frac{-b}{120}`	`frac{1}{6!}*e`
`c*f_i`	c	0	0	0	0	0	0
`d*f_(i+1)`	d	d	`frac{d}{2}`	`frac{d}{6}`	`frac{d}{24}`	`frac{d}{120}`	`frac{d}{6!}*e`
`e*f_(i+2)`	e	`2*e`	`2*e`	`frac{8*e}{6}`	`frac{16*e}{24}`	`frac{32*e}{120}`	`frac{2^6}{6!}*e`
	0	0	1	0	0	?	?

Considered five relations , since we required five unknown values using ,matlab

close all
clear all
clc
a = [1 1 1 1 1 ; -2 -1 0 1 2 ; 2 1/2 0 1/2 2 ; -4/3 -1/6 0 1/6 4/3 ; 2/3 1/24 0 1/24 2/3];
b = [0 ; 0 ; 1 ; 0 ; 0];
x = inv(a)*b;

Got values as ,

a = -0.0833
b =  1.3333
c = -2.5000
d =  1.3333
e = -0.0833

`frac{del^2f}{delx^2} approx a*f(x_(i-2))+b*f(x_(i-1))+c*f(x_i)+d*f(x_(i+1))+e*f(x_(i+2)) approx (-2*a -b + d + 2*d)*f^2(x)*h^2 + O(h^6)` --> (7)

On simple manipulating the equation by using all equation data about , we got equation (7)

`frac{del^2f}{delx^2} approx frac{-0.0833*f(x_(i-2))+1.3333*f(x_(i-1))+-2.5000*f(x_i)+1.3333*f(x_(i+1))+0.0833*f(x_(i+2))}{(Deltax)^2} approx (-2*0.0833 + 1/2*1.3333 + 1/2*1.3333 + 2*0.0833)*f^2(x) + O(h^4)``

`frac{del^2f}{delx^2} approx frac{-0.0833*f(x_(i-2))+1.3333*f(x_(i-1))+-2.5000*f(x_i)+1.3333*f(x_(i+1))+0.0833*f(x_(i+2))}{(Deltax)^2} approx (1.3333)*f^2(x) + O(h^4)``

Forward Difference (or) Skewed right scheme :

Same approach for any scheme , either central or skewed.

Here , we take 6 nodes(4+2) , all aligned in right side.

`frac{del^2f}{delx^2} approx a*f(i) + b*f(i+1) + c*f(i+2) + d*f(i+3) + d*f(i+4) + e*f(i+5)` ----> (2a)

close all
clear all
clc
a = [1 1 1 1 1 1 ; 0 1 2 3 4 5  ; 0 1/2 2 9/2 8 25/2  ; 0 1/6 4/3 9/2 32/3 125/6 ; 0 1/24 2/3 27/8 32/3 625/24 ; 0 1/120 32/120 243/120 1024/120 3125/120];
b = [0 ; 0 ; 1 ; 0 ; 0 ; 0];
x = inv(a)*b;

Results :

a = 3.7500
b = -12.8333
c = 17.8333
d = -13.0000
e = 5.0833
e = -0.83333

Substitute values in equaiton (2a)

`frac{del^2f}{delx^2} approx frac{3.7500*f(i) - 12.8333*f(i+1) + 17.8333*f(i+2) -13.0000*f(i+3) + 5.0833*f(i+4) - 0.8333*f(i+5)}{(1.0001*dx^2)}`

Backward Difference (or) Skewed-Left scheme :

`frac{del^2f}{delx^2} approx a*f(i) + b*f(i-1) + c*f(i-2) + d*f(i-3) + d*f(i-4) + e*f(i-5)` ` -----> (3a)

%%
clc
a3 = [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 -2^3/6 -3^3/6 -4^3/6 -5^3/6; b30 1/24 2^4/24 3^4/24 4^4/24 5^4/24 ;0 -1/120 -2^5/120 -3^5/120 -4^5/120 -5^5/120 ];
c3 = [0; 0; 1; 0; 0; 0];
b3 = inv(a3)*c3;

Results :

a = 3.7500
b = -12.8333
c = 17.8333
d = -13.0000
e = 5.0833
e = -0.83333

`frac{del^2f}{delx^2} approx frac{3.7500*f(i) - 12.8333*f(i-1) + 17.8333*f(i-2) -13.0000*f(i-3) + 5.0833*f(i-4) - 0.8333*f(i-5)}{(1.0001*dx^2)}` `

Main code snippet:

%%
close all
clear all
clc
% central , skewed-left , skewed right appproimations of second order
% derivative
% considered target node , x = pi/3;
x = pi/3;
% Grid values
dx = linspace(pi/4,pi/4000,20);
for i = 1:length(dx)
           [sk_ri(i)] = skewed_right(x,dx(i));
           [sk_le(i)] = skewed_left(x,dx(i));
           [ce_Ap(i)] = central_Approx_Taylor(x,dx(i));
end
loglog(dx,sk_ri)
hold on
loglog(dx,sk_le)
loglog(dx,ce_Ap)
xlabel('grid sizen on loglog scale')
ylabel('error values')
legend('skewed-right','skewed-left','central-approximation')
grid on

figure(2)
plot(log10(dx),log10(sk_ri),'color','b')
hold on
plot(log10(dx),log10(sk_le),'color','r')
plot(log10(dx),log10(ce_Ap),'color','m')
xlabel('grid sizen on log10 scale')
ylabel('error values')
legend('skewed-right','skewed-left','central-approximation')
grid on

figure(3)
plot(dx,sk_ri,'color','r')
plot(dx,sk_le,'color','b')
plot(dx,ce_Ap,'color','m')
xlabel('grid size')
ylabel('error values')
legend('skewed-right','skewed-left','central-approximation')
grid on

Central_Approx_code_Snippet:

function [out_central] =central_Approx_Taylor(x,dx)

out_fourth = ((-0.0833*dis_r(x-2*dx) + 1.3333*dis_r(x-dx)-2.5000*dis_r(x) + 1.3333*dis_r(x+dx) + 0.0833-dis_r(x+2*dx)))/(1.3333*dx^2);
True_value =  -2*exp(x)*sin(x);
out_central = abs(True_value-out_fourth);

end

Skewed_left_Approx_code_Snippet:

function [out_left] = skewed_left(x,dx)

out_left1 = ((3.7500*dis_r(x) - 12.8333*dis_r(x-dx)+17.8333*dis_r(x-2*dx) -13.0000*dis_r(x-3*dx) + 5.0833*dis_r(x-4*dx) - 0.8333*dis_r(x-5*dx)))/(1.0001*dx^2);
True_value = -2*exp(x)*sin(x);
out_left = abs(True_value-out_left1);

end

Skewed_Right_Approx_code_Snippet:

function [out_right] = skewed_right(x,dx)

out_fourth = ((3.7500*dis_r(x) - 12.8333*dis_r(x+dx)+17.8333*dis_r(x+2*dx) -13.0000*dis_r(x+3*dx) + 5.0833*dis_r(x+4*dx) - 0.8333*dis_r(x+5*dx)))/(1.0001*dx^2);
True_value = -2*exp(x)*sin(x);
out_right = abs(True_value-out_fourth);

end

Function :

function [z] = dis_r(x)
  
    z = exp(x)*cos(x);
end

Explanation of code :

1. Assigned required variables for computations -> Target node and grid size.

2. Declared Functions for all 3 schemes for second order derivative.

3. Using these function got the error for all schemes for different grid sizes.

4. plotted the data using loglog command.

As observed from the graph , central Approximation giving more error , for second order derivative for fourth order approximation , may be due to round-off erros in computations.

Ceantral Approximation are not useful in some applications , when the path of study in unidirectional , like Heat tranfer from a boundary.

In this case we can't use central Approximation.

So , we should consider the numerical scheme carefully , that depends on approximation order and h reinfinement.