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

CHALLENGE REPORT:

AIM:

To 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 ( Central difference, Skewed right-sided difference, Skewed left-sided difference).

OBJECTIVE:

* To derive the following 4th order approximations of the second-order derivative for Central difference, Skewed right-sided difference, Skewed left-sided difference.

* 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.

* Provide a plot that compares the absolute error between the above-mentioned schemes.

THEORY :

SCHEME DERIVATIONS:

CENTRAL DIFFERENCING SCHEME:

NO OF NODES:

The number of nodes = p+q-1

Where,

p is order of derivative.

q is order of approximations.

The no of nodes = 2+4-1=5.

DERIVATION:

To derive a fourth order approximation for the second order derivative,we start with these approximation,

`(∂ ^2 u)/ (∂ x^ 2) = (a f ( i − 2 ) + b f ( i − 1 ) + c f ( i ) + d f ( i + 1 ) + e f ( i + 2 ))/( d x^2)`

taylor series derivation for the above terms,

`af(i−2)=af(i)−(2af'(i)dx)/1!+(4af''(i)dx^2)/2!-(8af'''(i)dx^3)/3!+(16af''''(i)dx^4)/4!`

Simillarly expand all the terms in the equation (1) and are represented in the form of table using taylor table method.

`d x^ 2 f ' ' ( i )`	f(i)	dxf'(x)	`dx^2f''(i)`	`dx^3f'''(i)`	`dx^4f''''(i)`
af(i-2)	a	-2a	4a/2	-8a/6	16a/24
bf(i-1)	b	-b	b/2	-b/6	b/24
cf(i)	c	0	0	0	0
df(i+1)	d	d	d/2	d/6	d/24
ef(i+2)	e	2e	4e/2	8e/6	16e/24

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

A.X = B

X = inv(A)*B

we get the coefficients as,

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

substituting these coefficient in equation(1) we get,

`(∂ ^2 u)/ (∂ x ^2) = (( − 0.0833 ) f ( i − 2 ) + ( 1.3333 ) f ( i − 1 ) + ( − 2.5000 ) f ( i ) + ( 1.3333 ) f ( i + 1 ) + ( − 0.0833 ) f ( i + 2 ) )/(d x^2)`

SKEWED RIGHT HAND DIFFERENCE SCHEME:

NO OF NODES:

The number of nodes = p+q

Where,

p is order of derivative.

q is order of approximations.

The no of nodes = 2+4 =6.

DERIVATION:

To derive a fourth order approximation for the second order derivative,we start with these approximation,

`(∂ ^2 u)/( ∂ x^ 2) = (a f ( i ) + b f ( i + 1 ) + c f ( i + 2 ) + d f ( i + 3 ) + e f ( i + 4 ) + g f ( i + 5 ) )/(d x^2)`

expand all the terms in the equation (2) using taylor series expansion and are represented in the form of table using taylor table method.

`dx^2f''(x)`	f(i)	dxf'(i)	`dx^2f''(i)`	`dx^3f'''(i)`	`dx^4f''''(i)`	`dx^5f'''''(i)`
af(i)	a	0	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24	b/120
cf(i+2)	c	2c	4c/2	8c/6	16c/24	32c/120
df(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120
ef(i+4)	e	4e	16e/2	64e/6	256e/24	1024e/120
gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120
sum of co_efficients	0	0	1	0	0	0

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

A.X = B

X = inv(A)*B

we get the coefficients as,

a=3.7500, b=-12.8333, c=17.8333, d=-13.0000 ,e=5.0833, g=-0.8333

substituting these coefficient in equation(2) we get,

`(∂^ 2u)/( ∂ x^ 2) = (( 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 ) )/(d x^2)`

SKEWED LEFT HAND DIFFERENCE SCHEME:

NO OF NODES:

The number of nodes = p+q

Where, p is order of derivative.

q is order of approximations.

The no of nodes = 2+4 =6.

DERIVATION:

To derive a fourth order approximation for the second order derivative,we start with these approximation,

`(∂ ^2u)/(∂ x^2) =( a f ( i − 5 ) + b f ( i − 4 ) + c f ( i − 3 ) + d f ( i − 2 ) + e f ( i − 1 ) + g f ( i ))/( dx^2)`

expand all the terms in the equation (3) using taylor series expansion and are represented in the form of table using taylor table method.

`dx^2f''(i)`	f(i)	dxf'(i)	`dx^2 f''(i)`	`dx^3f'''(i)`	`dx^4f''''(i)`	`dx^5f'''''(i)`
af(i-5)	a	-5a	25a/2	−125a/6	625a/24	-3125a/120
bf(i-4)	b	-4b	16b/2	−64b/6	256b/24	-1024b/120
cf(i-3)	c	-3c	9c/2	−27c/6	81c/24	-243c/120
df(i-2)	d	-2d	4d/2	−8d/6	16d/24	-32d/120
ef(i-1)	e	-e	e/2	−e/6	e/24	-e/120
gf(i)	g	0	0	0	0	0
sum of co_efficients	0	0	1	0	0	0

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

A.X = B

X = inv(A)*B

we get the coefficients as,

a=-0.8333, b=5.0833, c=-13.0000, d=17.8333, e=-12.8333, g=3.7500

substituting these coefficient in equation(3) we get

`(∂ ^2 u)/( ∂ x^ 2) =( ( − 0.8333 ) f ( i − 5 ) + ( 5.0833 ) f ( i − 4 ) + ( − 13.0000 ) f ( i − 3 ) + ( 17.8333 ) f ( i − 2 ) + ( − 12.8333 ) f ( i − 1 ) + ( 3.7500 ) f ( i ))/( d x^ 2)`

WHY A SKEWED SCHEME IS USEFUL:

The system which doesn’t contain sufficient information on both side of the point of interest. In such cases we cannot use central differencing scheme so skewed scheme is useful in such cases.

WHAT CAN A SKEWED SCHEME DO THAT A CD SCHEME CANNOT DO:

CDS cannot not be employed at the boundary nodes due to unavailability of right and left nodes at the same time. In such cases the skewed differencing scheme is used.

PROGRAMES IN MATLAB:

clear all
close all
clc
% 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=inv(A)*B


% 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=inv(A)*B


% A=[1 1 1 1 1 1;-5 -4 -3 -2 -1 0;25/2 16/2 9/2 4/2 1/2 0;-125/6 -64/6 -27/6 -8/6 -1/6 0;625/24 256/24 81/24 16/24 1/24 0;-3125/120 -1024/120 -243/120 -32/120 -1/120 0]
% B=[0;0;1;0;0;0]
% X=inv(A)*B

x=pi/3;
dx= linspace(pi/4000,pi/40,20);
for i=1:length(dx)
    central_difference_error_forth(i)=forth_order_central_difference(x,dx(i))
    skewed_right_difference_error_forth(i)=forth_order_skewed_right_difference(x,dx(i))
    skewed_left_side_difference_error_forth(i)=forth_order_skewed_left_side_difference(x,dx(i))
end

%plot
figure(1)
loglog(dx,central_difference_error_forth,'color', 'r','linewidth',2)
hold on
loglog(dx,skewed_right_difference_error_forth,'color', 'b','linewidth',2)
hold on
loglog(dx,skewed_left_side_difference_error_forth,'color', 'g','linewidth',2)
xlabel('dx')
ylabel('forth order error')
title('dx vs error')
legend('central difference error forth','skewed right difference error forth','skewed left side difference error forth')

function central_difference_error=forth_order_central_difference(x,dx)
analytical_derivative=(-2)*exp(x)*sin(x);
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=inv(A)*B;
central_difference_forth_order=((X(1)*exp(x-(2*dx))*cos(x-(2*dx)))+(X(2)*exp(x-dx)*cos(x-dx))+(X(3)*exp(x)*cos(x))+(X(4)*exp(x+dx)*cos(x+dx))+(X(5)*exp(x+(2*dx))*cos(x+(2*dx))))/(dx)^2
central_difference_error=abs(central_difference_forth_order-analytical_derivative)
end

function skewed_right_difference_error=forth_order_skewed_right_difference(x,dx)
analytical_derivative=(-2)*exp(x)*sin(x);
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=inv(A)*B;
skewed_right_difference_forth_order=((X(1)*exp(x)*cos(x))+(X(2)*exp(x+dx)*cos(x+dx))+(X(3)*exp(x+(2*dx))*cos(x+(2*dx)))+(X(4)*exp(x+(3*dx))*cos(x+(3*dx)))+(X(5)*exp(x+(4*dx))*cos(x+(4*dx)))+(X(6)*exp(x+(5*dx))*cos(x+(5*dx))))/(dx)^2
skewed_right_difference_error=abs(skewed_right_difference_forth_order-analytical_derivative)
end

function skewed_left_side_difference_error=forth_order_skewed_left_side_difference(x,dx)
analytical_derivative=(-2)*exp(x)*sin(x);
A=[1 1 1 1 1 1;-5 -4 -3 -2 -1 0;25/2 16/2 9/2 4/2 1/2 0;-125/6 -64/6 -27/6 -8/6 -1/6 0;625/24 256/24 81/24 16/24 1/24 0;-3125/120 -1024/120 -243/120 -32/120 -1/120 0];
B=[0;0;1;0;0;0];
X=inv(A)*B;
skewed_left_side_difference_forth_order=((X(1)*exp(x-(5*dx))*cos(x-(5*dx)))+(X(2)*exp(x-(4*dx))*cos(x-(4*dx)))+(X(3)*exp(x-(3*dx))*cos(x-(3*dx)))+(X(4)*exp(x-(2*dx))*cos(x-(2*dx)))+(X(5)*exp(x-dx)*cos(x-dx))+(X(6)*exp(x)*cos(x)))/(dx)^2
skewed_left_side_difference_error=abs(skewed_left_side_difference_forth_order-analytical_derivative)
end

CODE EXPLANATION:

* clear all, close all, clc is used to clear all the variables and delet the graph, plot which are created already.

* Input value is assigned to variables x and range operator is used to assign dx value (dx range = [pi/4000,pi/40]).

* Use for loop with i ranges from 1 to length of dx and assign the central and skewed right differance and skewed left differance function which is created.

* A function is created for 4th order error with skewed right side difference error as the function definition and forth order skewed right side difference as the function name and x,dx as the input and assign the skewed right side difference approximation and analytical differentiate formulas and end the function.

* Create next function for 4th order error with central difference error as the function definition and forth order central differance as the function name and x,dx as the input and assign the central difference approximation and analytical differentiate formulas and end the function.

* Then create another function for 4th order error with skewed left side difference error as the function definition and forth order skewed left side differenceas the function name and x,dx as the input and assign the skewed left side difference approximation and analytical differentiate formulas and end the function.

* To plot the graph loglog command is used for (dx vs central difference error , dx vs skewed right side difference error and dx vs skewed left side difference error) which is used to compare the curve in detaily and xlabel and ylabel is added and legend is also added.

CONCLUSION:

* From the graph we can conclude that the central differencing scheme has very less error compared to other schemes. Therefore, central differencing is better than other schemes because it uses information from both the sides from the point of consideration to give approximative and accurate values.

* As the left and right skewed difference methods uses data only from one side of point of interest it gives less approximative and accurate values about the point.