Taylor Table Method using Matlab Code

Aim : To derive fourth-order approximation equation of second-order derivative using Matlab & compare it with Analytical Method.

Methods used:

Analytical:

Given expression :

f(x) = e^x*cos(x)

f’(x) = e^x(cos(x)-sin(x))

f”(x) = -2*e^xsin(x) (second derivative of the equation)

Numerical:`f^'''(i)*∆x^3`

Central Differencing :

The number of terms required is determined by this formula:

n = p+q-1

where p = order of approximation = 4

q = order of derivative = 2

so here, n = 5

Using Taylor’s Series,

(d^2 u)/〖dx〗^2 *∆x^2=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)

a*f(1-2) = a*f(i)- 〖a*f'(i)*(2∆x)〗^1/1! + 〖a*f''(i)*(2∆x)〗^2/2! - 〖a*f'''(i)*(2∆x)〗^3/3! + 〖a*f''''(i)*(2∆x)〗^4/4! - 〖a*f'''''(i)*(2∆x)〗^5/5! + - - - -

Applying Taylor Table Method:

	f(i)	f^' (i)*∆x	f^''(i)*∆x^2	f^'''(i)*∆x^3	f^''''(i)*∆x^4
a*f(i-2)	a	-2*a	4*a/2	-8*a/6	16*a/24
b*f(i-1)	b	-b	b/2	-b/6	b/24
c*f(i)	c	0	0	0	0
d*f(i+1)	d	d	d/2	b/6	d/24
e*f(i+2)	e	2*e	4*e/2	8*e/6	16*e/24

Algebraic Equation :

a+b+c+d+e = 0

-2a-b+d+2e = 0

4a/2+b/2+d/2+4e/2 = 1

-8a/6-b/6+b/6+8e/6 = 0

16a/24+b/24+d/24+16e/24 = 0

Evaluating by matrix inversion method in matlab X = A\B

a = -0.0833

b = 1.3333

c = 2.5000

d = 1.3333

e = -0.0833

Central differencing function code:

function output = Central_difference(x,dx)

%Analytical Derivative

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

%Central 

F = [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];

G = [0; 0; 1; 0; 0]

Z = F\G

a = Z(1);
b = Z(2);
c = Z(3);
d = Z(4);
e = Z(5);

%Numerical Derivative

a_into_i_minus_2 = a*(exp(x-(2*dx))*cos(x-(2*dx)));
b_into_i_minus_1 = b*((exp(x-dx))*cos(x-dx));
c_into_i = c*(exp(x)*cos(x));
d_into_i_plus_1 = d*((exp(x+dx))*(cos(x+dx)));
e_into_i_plus_2 = e*(exp(x+(2*dx))*cos(x+(2*dx)));

N_D = ((a_into_i_minus_2)+(b_into_i_minus_1)+(c_into_i)+(d_into_i_plus_1)+(e_into_i_plus_2))/(dx^2);

output = abs(A_D-N_D);

end

Skewed Left Sided :

The number of terms required is determined by this formula:

n = p+q

where p = order of approximation = 4

q = order of derivative = 2

so here, n = 6

Using Taylor’s Series,

(d^2 u)/〖dx〗^2 *∆x^2=a*f(i-1)+b*f(i-2)+c*f(i-3)+d*f(i-4)+e*f(i-5)+g*f(i-6)

Applying Taylor Table Method:

	f(i)	f^' (i)*∆x	f^''(i)*∆x^2	f^'''(i)*∆x^3	f^''''(i)*∆x^4	f^'''''(i)*∆x^5
a*f(i-1)	a	0	0	0	0	0
b*f(i-2)	b	-b	b/2	-b/6	b/24	-b/120
c*f(i-3)	c	-2*c	4*c/2	-8*c/6	16*c/24	-32*c/120
d*f(i-4)	d	-3*d	9*d/2	-27*d/6	81*d/24	-243*d/120
e*f(i-5)	e	-4*e	16*e/2	-64*e/6	256*e/24	-1024*e/120
g*f(i-6)	g	-5*g	25*g/2	-125*g/6	625*g/24	-3125*g/120

Algebraic Equation :

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

Evaluating by matrix inversion method in matlab X = A\B

a = 3.7500

b = -12.8333

c = 17.8333

d = -13.0000

e = 5.0833

g = -0.8333

Code to Skewed Left sided function:

function error_Left = skewed_lefted(x,dx)

%Analytical Derivative

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

%Skewed Left Sided

D = [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]

E = [0; 0; 1; 0; 0; 0]

Y = D\E

y1 = Y(1);
y2 = Y(2);
y3 = Y(3);
y4 = Y(4);
y5 = Y(5);
y6 = Y(6);

%Numerical Derivative

a_into_i = y1*(exp(x)*cos(x));
b_into_i_minus_1 = y2*((exp(x-dx))*cos(x-dx));
c_into_i_minus_2 = y3*(exp(x-(2*dx))*cos(x-(2*dx)));
d_into_i_minus_3 = y4*(exp(x-(3*dx))*cos(x-(3*dx)));
e_into_i_minus_4 = y5*(exp(x-(4*dx))*cos(x-(4*dx)));
f_into_i_minus_5 = y6*(exp(x-(5*dx))*cos(x-(5*dx)));

N_D = ((a_into_i)+(b_into_i_minus_1)+(c_into_i_minus_2)+(d_into_i_minus_3)+(e_into_i_minus_4)+(f_into_i_minus_5))/(dx^2);

error_Left = abs(A_D-N_D);

end

Skewed Right Sided :

The number of terms required is determined by this formula:

n = p+q

where p = order of approximation = 4

q = order of derivative = 2

so here, n = 6

Using Taylor’s Series,

(d^2 u)/〖dx〗^2 *∆x^2=a*f(i+1)+b*f(i+2)+c*f(i+3)+d*f(i+4)+e*f(i+5)+g*f(i+6)

Applying Taylor Table Method:

	f(i)	f^' (i)*∆x	f^''(i)*∆x^2	f^'''(i)*∆x^3	f^''''(i)*∆x^4	f^'''''(i)*∆x^5
a*f(i+1)	a	0	0	0	0	0
b*f(i+2)	b	b	b/2	b/6	b/24	b/120
c*f(i+3)	c	2*c	4*c/2	8*c/6	16*c/24	32*c/120
d*f(i+4)	d	3*d	9*d/2	27*d/6	81*d/24	243*d/120
e*f(i+5)	e	4*e	16*e/2	64*e/6	256*e/24	1024*e/120
g*f(i+6)	g	5*g	25*g/2	125*g/6	625*g/24	3125*g/120

Algebraic Equation :

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

Evaluating by matrix inversion method in matlab X = A\B

a = 3.7500

b = -12.8333

c = 17.8333

d = -13.0000

e = 5.0833

g = -0.8333

Code to skewed Right sided function:

function error_right = skewed_righted(x,dx)

%Analytical Derivative

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

%Skewed Right Sided

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 = A\B

x1 = X(1);
x2 = X(2);
x3 = X(3);
x4 = X(4);
x5 = X(5);
x6 = X(6);

%Numerical Derivative

a_into_i = x1*(exp(x)*cos(x));
b_into_i_plus_1 = x2*((exp(x+dx))*cos(x+dx));
c_into_i_plus_2 = x3*(exp(x+(2*dx))*cos(x+(2*dx)));
d_into_i_plus_3 = x4*(exp(x+(3*dx))*cos(x+(3*dx)));
e_into_i_plus_4 = x5*(exp(x+(4*dx))*cos(x+(4*dx)));
f_into_i_plus_5 = x6*(exp(x+(5*dx))*cos(x+(5*dx)));

N_D = ((a_into_i)+(b_into_i_plus_1)+(c_into_i_plus_2)+(d_into_i_plus_3)+(e_into_i_plus_4)+(f_into_i_plus_5))/(dx^2);

error_right = abs(A_D-N_D);

end

We know that all the coefficients of f’(i), f’’’(i), f’’’’(i), f’’’’’(i) are 0 except for f’’(i) which is 1.

Dividing both sides by ,

the truncation error is of order 4 thereby making this the 4^th order approximate

Main program:

clear all 
close all
clc

x = pi/3

dx = linspace(pi/10, pi/4000,100);

for i = 1:length(dx)

   Central_difference_error(i) = Central_difference(x,dx(i));
   Skewed_right_sided_error(i) = skewed_righted(x,dx(i));
   Skewed_left_sided_error(i) = skewed_lefted(x,dx(i));
   
end

loglog(Central_difference_error,'-*','color','r','markersize',10,'linewidth',2)
hold on
loglog(Skewed_right_sided_error,'-O','color','b','markersize',10,'linewidth',2)
hold on
loglog(Skewed_left_sided_error,'-<','color','g','markersize',10,'linewidth',2)
ylabel('error')
legend('Central_difference_error','Skewed_right_sided_error','Skewed_left_sided_error')

Plot:

From the plot, the error in CDS is lower compared to the skewed left & the right sided ones. The error decreases with decrease in grid size upto certain value & it increases again making it the best for the better approximate result.

why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

Skewed scheme take their data only from one direction whatever be the number of points available on the grid. For calculating corner node with higher approximation we need this method unlike CDS which takes it’s data from either sided.

Ex : Heat flow in a pipe, heat flows only in one direction, skewed schemes are preferred for faster, accurate flow rate derivatives.

Conclusion:

According to the flow, the three schemes should be applied where & when required with proper conditions to obtain the accurate results.