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

AIM:

To derive the fourth order approximations of a second order derivative using central differencing scheme, skewed right sided difference and skewed left sided difference with the help of taylor table method and to compare the analytical derivative with obtained values for the given function.

OBJECTIVE

• 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 that have derived.
• Provide a plot that compares the absolute error between the above-mentioned schemes
• In your own words, describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?

THEORY

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,

`(∂^2u)/(∂x2)=(af(i−2)+bf(i−1)+cf(i)+df(i+1)+ef(i+2))/(dx^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.

`dx^2f''(i)`	f(i)	`f'(i)dx`	`f''(i)dx2`	`f'''(i)dx3`	`f''''(i)dx4`
af(i-2)	a	-2a	2a	`-(4a)/3`	`-(2a)/3`
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	2e	`(4e)/3`	`(2e)/3`
sum of co_efficients	0	0	1	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.0833, b=1.3333, c=-2.5000, d=1.3333 e=-0.0833

substituting these coefficient in equation(1) we get,

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

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,

`(∂^2u)/(∂x^2)=(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/(dx^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''(i)`	f(i)	`f'(i)dx`	`f''(i)dx^2`	`f'''(i)dx^3`	`f''''(i)dx^4`	`f'''''(i)dx^5`
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))/(dx^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)=(af(i−5)+bf(i−4)+cf(i−3)+df(i−2)+ef(i−1)+gf(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)	`f'(i)dx`	`f''(i)dx^2`	`f'''(i)dx^3`	`f''''(i)dx^4`	`f'''''(i)dx^5`
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

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

STEPS

• First we created taylor's table for all three schemes
• Then we created linear eqautions
• The linear equation obtained from the above table are solved using matrix inversion method
• By using the command window we finds the values of coefficients for all three schemes
• Then we substitute those values in taylor series to do approximations
• In the main code we called functions (which are seperately created ) in for loop to do approximation for varying value of dx

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

We created three seperate function codes for three approximations

Function 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

Function code skewed right hand approximation

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

Function code skewed left hand approximation

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);

OUTPUT

OBSERVATIONS

• As the order of approximation increases error decreases
• For certain range the value of error for all three approximation is almost constant
• After some point the value of error changes
• Skewed left hand approximation have some diiferent graph
• Central approximation have higher value of error
• Skewed right has less value of error

CONCLUSION

WHY A SKEWED SCHEME IS USEFUL

The central differencing scheme has some limitation.If there no sufficient information on side of inteterst,we cannot adopt CDST .So skewed scheme is useful in such cases.

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

CDS can't 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

matlab files is attached below