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

AIM:- Derive The 4th order Approximation for 2nd order Derivative, and find the right scheme.

OBJECTIVE:-

• Derive the following 4th order approximations of the second-order derivative.

1. Central difference
2. Skewed right-sided difference
3. Skewed left-sided difference

• After deriving these schemes, prove that your skewed schemes are fourth-order accurate
• 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 you have derived.
• compares the absolute error between the above-mentioned schemes.
• describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?

PROCEDURE:-

1) Central Difference Scheme:

we are going to derive Central difference Scheme

let we have to calculate Number of nodes

N= P+Q-1

Where, N= No. of nodes

P= order of derivative

Q= order of Approximation.

Now,

N=2+4-1

N=5;

So Numerical Stencil for Central Difference is as follows

[i-2]---------[i-1]---------[i]---------[i+1]--------[i+2]

Second order derivative for central difference Scheme will be Written as

`(d^2y)/(dx^2)=af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`----------------eq (1)

By using taylors series Expansion we can Written above equation as,

`af(i-2)=a(f(i)-2Deltaxf'(i)+(2Deltax^2f''(i))/(2!)-(2Deltax^3f'''(i))/(3!)+(2Deltax^4f^''''(i))/(4!))`

`af(i-1)=b(f(i)-Deltaxf'(i)+(Deltax^2f''(i))/(2!)-(Deltax^3f'''(i))/(3!)+(Deltax^4f^''''(i))/(4!))`

`cf(i)=cf(i)`

`df(i+1)=d(f(i)+Deltaxf'(i)+(Deltax^2f''(i))/(2!)+(Deltax^3f'''(i))/(3!)+(Deltax^4f^''''(i))/(4!))`

`ef(i+2)=e(f(i)+2Deltaxf'(i)+(2Deltax^2f''(i))/(2!)+(2Deltax^3f'''(i))/(3!)+(2Deltax^4f^''''(i))/(4!))`

we can find Coeefficient value by Taylors table

	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`
af(i-2)	a	-2a	(4/2)*a	-(8/6)*a	(16/24)*a
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	(4/2)*e	(8/6)*e	(16/24)*e
Sum of Coeifficients	0	0	1	0	0

Table.1 Taylor table for Central Difference Scheme

Form Taylor table we got the equation

a+b+c+d+e=0

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

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

-(4/3)*a-b/6+d/6+(4/3)*e=0

(2/3)*a+b/24+d/24+(2/3)*e=0

converting above equation in matrix form

A1=[1 1 1 1 1;-2 -1 0 1 2;2 1/2 1/2 2;-4/3 1/6 1/6 4/3;2/3 1/24 1/24 2/3]

B1=[0;0;1;0;0]

X1=[a;b;c;d;e]

`therefore` A1X1=B1

by solving matrix inversion method in matlab we can get,

a=-0.0833

b=1.3333

c=-2.5000

d=1.3333

e=-0.0833

hence a=d and b=d;

equation 1 becomes,

`af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)=Deltax^2f''(i)`

`f''(i)=(af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/(Deltax^2)`

This is the expression for Central Differrence Scheme.

2) Skewed Right hand Scheme:

we are going to Skewed Right hand Scheme

let we have to calculate Number of nodes

N= P+Q

Where, N= No. of nodes

P= order of derivative

Q= order of Approximation.

Now,

N=2+4

N=6;

So Numerical Stencil for Skewed Right hand Scheme is as follows

[i]---------[i+1]---------[i+2]---------[i+3]--------[i+4]----------[i+5]

for second Derivative for Skewed Right hand Scheme can be written as

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

By using taylors series Expansion we can Written above equation as,

`af(i)=af(i)`

`bf(i+1)=b(f(i)+Deltaxf'(i)+((Deltax^2)f''(i))/(2!)+((Deltax^3)f'''(i))/(3!)+((Deltax^4)f''''(i))/(4!)+((Deltax^5)f'''''(i))/(5!))`

`cf(i+2)=c(f(i)+2Deltaxf'(i)+((2Deltax^2)f''(i))/(2!)+((2Deltax^3)f'''(i))/(3!)+((2Deltax^4)f''''(i))/(4!)+((2Deltax^5)f'''''(i))/(5!))`

`df(i+3)=d(f(i)+3Deltaxf'(i)+((3Deltax^2)f''(i))/(2!)+((3Deltax^3)f'''(i))/(3!)+((3Deltax^4)f''''(i))/(4!)+((3Deltax^5)f'''''(i))/(5!))`

`ef(i+4)=e(f(i)+4Deltaxf'(i)+((4Deltax^2)f''(i))/(2!)+((4Deltax^3)f'''(i))/(3!)+((4Deltax^4)f''''(i))/(4!)+((4Deltax^5)f'''''(i))/(5!))`

`gf(i+5)=g(f(i)+5Deltaxf'(i)+((5Deltax^2)f''(i))/(2!)+((5Deltax^3)f'''(i))/(3!)+((5Deltax^4)f''''(i))/(4!)+((5Deltax^5)f'''''(i))/(5!))`

we can find Coeefficient value by Taylors table

	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`	`Deltax^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	2c	(4/3)*c	(16/24)*c	(32/120)*c
df(i+3)	d	3d	(9/2)*d	(9/2)*d	(81/24)*d	(243/120)*d
ef(i+4)	e	4e	8e	(32/3)*e	(256/24)*e	(1021/120)*e
gf(i+5)	g	5g	(25/2)*g	(125/6)*g	(625/24)*g	(3125/120)*g
Sum of coefficient	0	0	1	0	0	0

Table.2 Taylor table for Skewed Right hand Scheme

From taylor table we can get,

a+b+c+d+e+g=0

b+2c+3d+4e+5g=0

(b/2)+2c+(9/2)d+8e+(25/2)g=1

(b/24)+(16/24)c+(81/24)d+(256/24)e+(625/24)g=0

(b/120)+(32/120)c+(243/120)d+(1021/120)e+(3125/120)g=0

Converting above equation in Matrix form,

A2=[1 1 1 1 1 1; 0 1 2 3 4 5; 0 1/2 2 9/2 8 25/2; 0 1/24 16/24 81/24 256/24 625/24; 0 1/120 32/120 243/120 1021/120 3125/120];

B2=[0;0;1;0;0;0];

X2=[a,b,c,d,e,g];

by considering A2X2=B2;

by solving above equation in Matlab by inversion method we get,

a=3.7500

b=-12.8333

c=17.8333

d=-13.0000

e=5.0833

g=-0.8333

we having,

af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)=`Deltax^2f''(i)`

f''(i)=`(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/(Deltax^2`

this is the expression for Skewed Right hand Scheme.

3) Skewed Left hand Scheme:

we are going to Skewed Left hand Scheme

let we have to calculate Number of nodes

N= P+Q

Where, N= No. of nodes

P= order of derivative

Q= order of Approximation.

Now,

N=2+4

N=6;

So Numerical Stencil for Skewed Right hand Scheme is as follows

[i-5]---------[i-4]---------[i-3]---------[i-2]--------[i-1]----------[i]

for second Derivative for Skewed Right hand Scheme can be written as

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

By using taylors series Expansion we can Written above equation as,

`af(i-5)=a(f(i)-5Deltaxf'(i)+((5Deltax^2)f''(i))/(2!)-((5Deltax^3)f'''(i))/(3!)+((5Deltax^4)f''''(i))/(4!)-((5Deltax^5)f'''''(i))/(5!))`

`bf(i-4)=b(f(i)-4Deltaxf'(i)+((4Deltax^2)f''(i))/(2!)-((4Deltax^3)f'''(i))/(3!)+((4Deltax^4)f''''(i))/(4!)-((4Deltax^5)f'''''(i))/(5!))`

`cf(i-3)=c(f(i)-3Deltaxf'(i)+((3Deltax^2)f''(i))/(2!)-((3Deltax^3)f'''(i))/(3!)+((3Deltax^4)f''''(i))/(4!)-((3Deltax^5)f'''''(i))/(5!))`

`df(i-2)=d(f(i)-Deltaxf'(i)+((Deltax^2)f''(i))/(2!)-((Deltax^3)f'''(i))/(3!)+((Deltax^4)f''''(i))/(4!)-((Deltax^5)f'''''(i))/(5!))`

`ef(i-1)=e(f(i)-Deltaxf'(i)+((Deltax^2)f''(i))/(2!)-((Deltax^3)f'''(i))/(3!)+((Deltax^4)f''''(i))/(4!)-((Deltax^5)f'''''(i))/(5!))`

`gf(i)=gf(i)`

we can find Coeefficient value by Taylors table

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

Table.3 Taylor table for Skewed Left hand Scheme

From taylor table we can get,

a+b+c+d+e+g=0

-5a-4b-3c-2d-e=0

(25/2)a+(16/2)b+(9/2)c+(4/2)d+(e/2)=1

(-125/6)a+(-64/6)b+(-27/6)c+(-8/6)d+(-e/6)=0

(625/24)a+(256/24)b+(81/24)c+(16/24)d+(e/24)=0

(-3125/120)a+(-1024/120)b+(-243/120)c+(-32/120)d+(-e/120)=0

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

B3=[0;0;1;0;0;0];

X3=[a,b,c,d,e,g];

by considering A3X3=B3;

by solving above equation in Matlab by inversion method we get,

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

we having,

`f''(i)=(af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i))/(Deltax^2)`=`Deltax^2f''(i)`

This is the expression for Skewed Left hand Scheme.

MATLAB CODE:-

clear all
close all
clc

% Analytical Function=exp(x)*cos(x)
x=pi/3;

% Second derivative of the function is 
second_derivative_Analytical=(-2)*exp(x)*sin(x);

% calculating the Coefficients of the central Difference Scheme

A1=[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];
B1=[0;0;1;0;0];
X1=linsolve(A1,B1);

% calculating the Coefficients of the Skewed Right hand Scheme

A2=[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];
B2=[0;0;1;0;0;0];
X2=linsolve(A2,B2);

% calculating the Coefficients of the Skewed Left hand Scheme

A3=[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];
B3=[0;0;1;0;0;0];
X3=linsolve(A3,B3);

% Numerical derivative for Given Schemes are

dx=linspace(pi/40,pi/4000,20);
for i=1:length(dx)
    % for Central Difference Scheme
     Central_Difference(i)=(((X1(1)*exp(x-(2*dx(i)))*cos(x-(2*dx(i)))))+(X1(2)*exp(x-(dx(i)))*cos(x-(dx(i))))+(X1(3)*exp(x)*cos(x))+(X1(4)*exp(x+(dx(i)))*cos(x+(dx(i))))+(X1(5)*exp(x+(2*dx(i)))*cos(x+(2*dx(i)))))/(dx(i)^2);

    % for  Skewed Right hand Scheme
    Right_Skewed(i)=(((X2(1)*exp(x)*cos(x))+(X2(2)*exp(x+dx(i)))*cos(x+dx(i))+((X2(3)*exp(x+(2*dx(i)))*cos(x+(2*dx(i)))))+(X2(4)*exp(x+(3*dx(i)))*cos(x+(3*dx(i))))+(X2(5)*exp(x+(4*dx(i))))*cos(x+(4*dx(i)))+(X2(6)*exp(x+(5*dx(i)))*cos(x+(5*dx(i))))))/(dx(i)^2);

    % for Skewed Left hand Scheme
    Left_Skewed(i)=((((X3(1)*exp(x-(5*dx(i))))*cos(x-(5*dx(i))))+(X3(2)*exp(x-(4*dx(i))))*cos(x-(4*dx(i)))+(X3(3)*exp(x-(3*dx(i))))*cos(x-(3*dx(i)))+(X3(4)*exp(x-(2*dx(i))))*cos(x-(2*dx(i)))+(X3(5)*exp(x-(dx(i))))*cos(x-(dx(i)))+(X3(6)*exp(x)*cos(x)))/(dx(i)^2));

    % calculate Absolute Error

    Error_Central_Difference(i)=abs(Central_Difference(i)-second_derivative_Analytical);
    Error_Right_Skewed(i)=abs(Right_Skewed(i)-second_derivative_Analytical);
    Error_Left_Skewed(i)=abs(Left_Skewed(i)-second_derivative_Analytical);

end
loglog(dx,Error_Central_Difference,'r')
hold on
loglog(dx,Error_Right_Skewed,'b')
hold on
loglog(dx,Error_Left_Skewed,'g')
hold on
legend('central difference scheme','Right hand Skewed Scheme','Left hand Skewed Scheme')
xlabel('dx')
ylabel('absolute Error')

Explaination of Code:-

Step-1:- First we have to calculate the Second order derivative of the given Function.

Step-2: Then we can calculate the coeifficient by using Taylor table Expansion for each scheme Central differencing, Right Skewed and Left Skewed Schemes by using matrix inverse interpolation Method, like X1,X2,and X3 for Central differencing, Right Skewed and Left Skewed Schemes Respectively.

Step-3: by using for loop we have to write code for Calculation of Numerical Derivative for Central differencing, Right Skewed and Left Skewed Schemes Respectively, our 'x' value will be [pi/3] and 'dx' value will be [pi/40,pi/4000,20].

Step-4:- we have to calculate the Absolute Error for Central differencing, Right Skewed and Left Skewed Schemes Respectively.

Step-5:- Then we have to plot graph by using of 'loglog' command for 'absolute error' vs 'dx' for compare the Central differencing, Right Skewed and Left Skewed Schemes.

RESULT:-

Plot of Relation of Absolute Error and dx

Fig.1 Plot of Relation of Absolute Error and dx

Above graph shows that how the behavior of Absolute Error of Differrent Schemes( Central diffrencing, Right Skewed and Left Skewed changes with respect to dx value.

CONCLUSION:

• The above graph shows that In central Difference scheme having less absolute error as compared to Right Skewed and Left Skewed.
• When order of approximation Increases the error will be decreases.
• The rate of absolute error falles is almost same in Right Skewed and Left Skewed Schemes but in central diffrence is diffrent and comparatively low.

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

• When we have to get the information from both the side we required Central diffrence Scheme. but for Time constraint and Boundary condition problem we can't Use CD.
• when we required information from Left or Right side of given reference point Skewed Schemes are used. we can get the information for boundary level problems by using Skewed Schemes.
• we can conclude that we can get less error from Central diffrence Schemes but we can not use CD at less information we have.