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

Aim : To solve the 4th order approxiamtion of 2nd order derivative using :

1. Central difference

2. Skewed right-sided difference

3. Skewed left-sided difference

1. Central difference method: We take information from both sides of a point.

	f(i)	f'(i)*`Deltax`	f''(i)* `Deltax^2`	f'''(i)* `Deltax^3`	f''''(i)*`Deltax^4`
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
	0	0	1	0	0

We create matrix for the above table to solve for a,b,c,d and e

A = `[[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]]`

B = `[[0],[0],[1],[0],[0]]`

X = rats(linsolve(A,B))

We use rats to convert decimal to fractions we get:

where: a = -1/12

b = 4/3

c = -5/2

d = 4/3

e = -1/12

2.Right skew method: Here we take values right side of the point.

	f(i)	f'(i)*`Deltax`	f''(i)* `Deltax^2`	f'''(i)* `Deltax^3`	f''''(i)*`Deltax^4`	f'''''(i)*`Deltax^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
	0	0	1	0	0	0

We create matrix for the above table to solve for a,b,c,d,e and g

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 = rats(linsolve(A,B))

We use rats to convert decimal to fractions we get:

3. Left skew method : Here we take points on the left side of the point.

	f(i)	f'(i)*`Deltax`	f''(i)* `Deltax^2`	f'''(i)* `Deltax^3`	f''''(i)*`Deltax^4`	f'''''(i)*`Deltax^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
	0	0	1	0	0	0

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 = rats(linsolve(A,B))

We use rats to convert decimal to fractions we get:

Steps in codeing :

1. Create a function for `f(x)=cosx*e(x)`

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

2. Insert the coefficients in a matrix and solve them using linsolve

% coefficients for Left skew method
clear all
close all
clc

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 = rats(linsolve(A,B))


%Coefficients for Right skew method

clear all
close all
clc

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 = rats(linsolve(A,B))

%coefficienst for Central difference method

clear all
close all
clc

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

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

X = rats(linsolve(A,B))

3. Create the functions for Central difference,Right skew and Left skew and substitute the values obtained from above program in it.

%central difference
function out = centraldiff(x,dx)
a=-1
b=16
c=-30
d=16
e=-1
analyticalsol= -2*exp(x)*sin(x);
central = (a*func(x-2*dx)+b*func(x-dx)+c*func(x)+d*func(x+dx)+e*func(x+2*dx))/(12*(dx^2));
out = abs(central-analyticalsol);
end

%Right skew
function out = rightskew(x,dx)

analyticalsol = -2*exp(x)*sin(x);
rightskew = (45*func(x)-154*func(x+dx)+214*func(x+2*dx)-156*func(x+3*dx)+61*func(x+4*dx)-10*func(x+5*dx))/(12*(dx^2))
out = abs(analyticalsol-rightskew);
end

%Left skew
function out = leftskew(x,dx)

analyticalsol = -2*exp(x)*sin(x);
leftskew= (45*func(x)-154*func(x-dx)+214*func(x-2*dx)-156*func(x-3*dx)+61*func(x-4*dx)-10*func(x-5*dx))/(12*(dx^2))
out = abs(analyticalsol-leftskew);
end

4. Now write main code involving all these functions and plot the graph.

Main code:

clear all
close all
clc

x = pi/3;
dx = linspace(pi/4,pi/400,100);


%central
for i = 1:length(dx)
    central_error(i) = centraldiff(x,dx(i));
end

%right
for i = 1:length(dx)
    right_error(i)=rightskew(x,dx(i));
end

%left
for i = 1:length(dx)
    left_error(i)=leftskew(x,dx(i));
end

%plotting
figure(1)
loglog(dx,central_error,'LineWidth',1,'Marker','pentagram','MarkerSize',2)
hold on
loglog(dx,right_error,'LineWidth',1,'Marker','pentagram','MarkerSize',2)
loglog(dx,left_error,'LineWidth',1,'Marker','pentagram','MarkerSize',2)
legend('central difference','right-skew difference','left-skew difference')
xlabel('dx')
ylabel('Error')
grid on

Plots and Results :

Conclusion :

1.Central difference method has least error among right skew and left skew method.

2. Left skew method error decreases as value of dx increases

3. As dx values decreases Right skew and left skew method both converge and has same errror.