Week 3 - Taylor table method and Matlab code

Introduction

This project is to derive an equation to calculate the 2nd order numerical derivative for central difference,skewed right sided and skewed left sided schemes of 4th order accuracy.

Also to find the 2nd order derivatives of exp(x)*cos(x) function and to compute the erro between its analytical and numerical derivative for all the 3 schemes

Theory

The forth order approximation for a second order derivative at a point is found out using 3 different methods.

1. Skewed right sided difference

2. Skewed left sided difference

3. Central difference

Skewed left sided difference

The derivative at a point is found using the information at points on the left hand side from the point on which we are calculating the derivative. The number of points or nodes required to calculate the derivative is found using the formula p+q for one sided difference and p+q-1 for central difference.

where p is the order of derivative(in this case 2).

q is the order of approximation(in this case 4).

Hence the total number of nodes required to calculate the derivative is 6 nodes.

Derivation

The algebraic equation used to solve the derivative of a second order function with an aproximation of 4th order is found using the Taylors Table. The partial derivative is equated to a hypthetical function as below.

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

The above terms are expanded using taylor series. For e.g

`b f(i-1)= b( f(i) -f'(i) trianglex+f''(i) 1^2/(2!)trianglex^2-f'''(i) 1^3/(3!)trianglex^3+f''''(i) 1^4/(4!)trianglex^4-f'''''(i) 1^5/(5!)trianglex^5)`

Then the taylor table is formed from the sum of all the terms expanded using taylor series.

	`f(i)`	`f'(i)trianglex`	`f '(i)trianglex^2`	`f'''(i)trianglex^3`	`f''''(i)trianglex^4`	`f'''''(i)trianglex^5`
`a f(i)`	a	0	0	0	0	0
`b f(i-1)`	b	-b	`1^2/(2!) b`	`-1^3/(3!) b`	`1^4/(4!) b`	`-1^5/(5!) b`
`c f(i-2)`	c	-2c	`2^2/(2!)c`	`-2^3/(3!)c`	`2^4/(4!)c`	`-2^5/(5!)c`
`d f (i-3)`	d	-3d	`3^2/(2!)d`	`-3^3/(3!)d`	`3^4/(4!)d`	`-3^5/(5!)d`
`e f (i-4)`	e	-4e	`4^2/(2!) e`	`-4^3/(3!) e`	`4^4/(4!) e`	`-4^5/(5!) e`
`g f(i-5)`	g	-5g	`5^2/(2!)g`	`-5^3/(3!)g`	`5^4/(4!)g`	`-5^5/(5!)g`
	0	0	1	0	0	0

Using MATLAB,the equations are written in matrix form

and the D matrix contains the values of the coefficients a,b,c,d,e,f,g.

The coefficients are substituted in the hypothetical functions and divided by delx^2 thereby providing a left sided skewed scheme.

`(del^2y)/(delx^2)=(3.7500 f(i-5dx)-12.8333 f(i-4dx)-17.8333f(i-3dx)-13.0000f(i-2dx)-5.0833(i-1dx)+-0.8333f(i))/(trianglex^2)`

Reason for 4th order accuracy

`(del^2y)/(delx^2)=f(i)[...]+f'(i)trianglex[...]+f '(i)trianglex^2[...] +f'''(i)trianglex^3[...]+f''''(i)trianglex^4[...]+f'''''(i)trianglex^5[...]+f''''''(i)trianglex^6[....]`

The 2nd order partial derivative is equated to all the sum of terms in the taylor table. We know that terms that contain functions of 1st,3rd,4th,5th order all equate to zero. To obtain the 2nd order differential equation the above entire equation is divided by `trianglex^2`and this reduces the 6th derivative of function to order `triangle x^4`.This is called the truncation error. Hence this scheme is called as 4th order accurate.

`(del^2y)/(delx^2)=f''(i)+O(triangle x^4)`

In the similar way the schemes for skewed right sided and central difference are found out

Skewed right sided difference

Here the numerical stencil consits of 5 points to the right of the point at which we want the 2nd order partial derivative.

`(del^2y)/(delx^2)=(3.7500 f(i)+12.833 f(i+dx)+17.833 f(i+2dx)+13 f(i+3dx)+5 f(i+4dx)+0.83f(i+5dx))/(triangledx^2)`

Central difference

Here the numerical stencil is symmetric taking information from 2 points before and after the point at which the derivative is calculated.

`(del^2y)/(delx^2)= (-0.0833f(i-2dx)+1.3333f(i-dx)-2.5 f(i)+1.333f(i+dx)-0.0833f(i+2dx))/(trianglex^2)`

Code

Taylor_Table.m

1.The main code that contains the analytical derivative and calculates the error between the analytical derivative and 3 schemes by calling the respective functions

1.Central_difference

2.left_sided

3.right_sided

2.The functions taken in 2 variables as input ,1. the point at which the derivative is calculated and 2. the grid spacing values.

3. A plot is plotted between dx and the truncation eror that is found at by the absolute difference between analytical and numerical derivative.

x=pi/3;
dx=linspace(pi/40,pi/4000,50);
analytical_derivative=-2*exp(x)*sin(x);

for i=1:length(dx)
Central_difference_error(i)=abs(analytical_derivative-central_difference(x,dx(i)));
Skewed_left_sided_error(i)=abs(analytical_derivative-left_sided(x,dx(i)));
Skewed_right_sided_error_(i)=abs(analytical_derivative-right_sided(x,dx(i)));
end
loglog(dx,Central_difference_error,dx,Skewed_left_sided_error,dx,Skewed_right_sided_error_);
legend('Central Difference','Skewed Left Sided','Skewed Right Sided');
xlabel('dx');
ylabel('Truncation Error');

func.m

This function contains the function for which the 2nd order derivative is calculated

function f1=func(x)
f1=exp(x)*cos(x);
end

Central_difference.m

This function contains the algebraic terms for calculating the 2nd order difference using a symmetric stencil.

function partial_derivative=central_difference(x,dx)
partial_derivative=(-0.0833*func(x-(2*dx))+1.3333*func(x-dx)-2.5*func(x)+1.3333*func(x+dx)-0.0833*func(x+(2*dx)))/(dx^2);
end

Left_sided.m

This function contains the algebraic terms for calculating the 2nd order difference using a left sided skeweds stencil.

function partial_derivative_2=left_sided(x,dx)
partial_derivative_2=(-0.8333*func(x-(5*dx))+5.0833*func(x-(4*dx))-13*func(x-(3*dx))+17.8333*func(x-(2*dx))-12.8333*func(x-dx)+3.75*func(x))/(dx^2);
end

Right_sided.m

This function contains the algebraic terms for calculating the 2nd order difference using a right sided skewed stencil.

function partial_derivative_3=right_sided(x,dx)
partial_derivative_3=((3.7500*func(x))-(12.8333*func(x+dx))+(17.8333*func(x+(2*dx)))-(13.0000*func(x+(3*dx)))+(5.0833*func(x+(4*dx)))-(0.8333*func(x+(5*dx))))/(dx^2);
end

MatrixInverse.m

This is used to find the co efficients of the corresponding schemes. The inverse of matrix is used to find the coefficient of matrix.

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

B=[0;0;1;0;0;0;];

A=Central_Difference\B;

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

C=Right_Skewed\B;

Left_Skewed=[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;];
D=Left_Skewed\B;
   

Plot

Trend

1. As the grid spacing becomes smaller and smaller, the truncation erro i.e the difference between the analytical and numerical derivative becomes smaller.

2. With increase in dx ,the error increases in magnitude with change in skewed right sided scheme being the highest.

3. The change in error for the central difference is very small ,with trend being constant.

4. In the left sided scheme,for increase in dx the error for some value decreases and then bgins to increase.

Why a skewed scheme is useful?

Say for example a fluid is flowing over a flat plate and you want to compute the velocity of fluid at a point on the plate. This happens to be the boundary. Now if we were to use the central differencing method to compute the velocity we would require values at the top of the boundary and at the bottom.

The velocity above the boundary can be found out, whereas the values of velocity below the boundary cannot be as it happens to be outside the domain of computation.

In cases like these we can use one sided difference like left or right skewed scheme to obtain values at the boundary.