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

USING THE CENTRAL DIFFERENCE METHOD:

IN central difference method numerical stencil will be

(i-2) (i-1) (i) (i+1) (i+2)

for fourth order approximation at point (i) we have to take information from both side right& left and the distance between them will be ∆x.

our hypothesis equation will be :

∂²u/∂x² = af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)...............(1)

to find the value of (a,b,c,d,e) we are gone a use taylor table series method.

(a)f(i-2)= (a)f(i)-(a)f'(i)(2∆x)/1! +(a)f''(i)(2∆x)^2/2! - (a)f'''(i)(2∆x)^3/3! +(a)f''''(i)(2∆x)^4/4! - (a)f'''''(i)(2∆x)^5/5! + (a)f''''''(i)(2∆x)^6/6!+.............

(b)f(i-1)= (b)f(i)-(b)f'(i)(∆x)/1! + (b)f''(i)(∆x)^2/2! - (b)f'''(i)(∆x)^3/3! + (b)f''''(i)(∆x)^4/4! - (b)f'''''(i)(∆x)^5/5! + (b)f''''''(i)(∆x)^6/6!+ ..............

(c)f(i)=(c)f(i)

(d)f(i+1)=(d)f(i)+(d)f'(i)(∆x)/1! - (d)f''(i)(∆x)^2/2! + (d)f'''(i)(∆x)^3/3! - (d)f''''(i)(∆x)^4/4! + (d)f'''''(i)(∆x)^5/5! - (d)f''''''(i)(∆x)^6/6!+..........

(e)f(i+2)= (e)f(i)+(e)f'(i)(2∆x)/1! - (e)f''(i)(2∆x)^2/2! + (e)f'''(i)(2∆x)^3/3! - (e)f''''(i)(2∆x)^4/4! + (e)f'''''(i)(2∆x)^5/5! - (e)f''''''(i)(2∆x)^6/6!+........

arranging above terms in taylor tyabel:

	f(i)	∆xf(i)	∆x^2f''(i)	∆xf^3'''(i)	∆x^4f''''(i)	∆x^5f'''''(i)	∆x^6f''''''(i)
(a)f(i-2)	a	-2a	4a/2	-8a/6	16a/24	-32a/120	64a/240
(b)f(i-1)	b	-b	b/2	-b/6	b/24	-b/120	b/240
(c)f(i)	c	0	0	0	0	0	0
(d)f(i+1)	d	d	d/2	d/6	d/24	d/120	d/240
(e)f(i+2)	e	2e	4e/2	8e/6	16e/24	-32e/120	64e/240
sum of above all	0	0	1	0	0	?	?

from above this we can write in matrix form AX=B

[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] [a b c d e ]=[0 0 1 0 0 ]

after solving above matrix the values we get

a=-0.083333

b=1.333333

c=-2.500000

d=1.333333

e=-0.83333

dividing both the side by equation (1) by (∆x^2) we get

∂²u/∂x² = af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2) / (∆x^2) %we can put the above vlues(a b c d e) in it

after this to get the second order deravative the coeficient of f''(i) has to be 1 so as to get 1 we have to ignore others or set it as 0.

here we can also see that a is equal to e and b is equal to d we can aslso deravie that and solve the equation and get to the final equation.

USING THE SKEWED RIGHT SIDED DIFFERENCE:

In this skewed right hand side difference method we have to define the five points on the right of (i).

the numerical stencil will be

(i) (i+1) (i+2) (i+3) (i+4) (i+5)

the hypothesis eqquation will be

∂²u/∂x² = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)

writing taylor series expansion for each term

af(i)=af(i)

(b)f(i+1)= (b)f(i)+(b)f'(i)(∆x)/1! + (b)f''(i)(∆x)^2/2! + (b)f'''(i)(∆x)^3/3! + (b)f''''(i)(∆x)^4/4! + (b)f'''''(i)(∆x)^5/5! + (b)f''''''(i)(∆x)^6/6!+.........

(c)f(i+2)=(c)f(i)+(c)f'(i)(2∆x)/1! + (c)f''(i)(2∆x)^2/2! + (c)f'''(i)(2∆x)^3/3! + (c)f''''(i)(2∆x)^4/4! + (c)f'''''(i)(2∆x)^5/5! + (c)f''''''(i)(2∆x)^6/6!+.......

(d)f(i+3)=(d)f(i)+(d)f'(i)(3∆x)/1! + (d)f''(i)(3∆x)^2/2! + (d)f'''(i)(3∆x)^3/3! + (d)f''''(i)(3∆x)^4/4! + (d)f'''''(i)(3∆x)^5/5! + (d)f''''''(i)(3∆x)^6/6!+..........

(e)f(i+4)= (e)f(i)+(e)f'(i)(4∆x)/1! + (e)f''(i)(4∆x)^2/2! + (e)f'''(i)(4∆x)^3/3! + (e)f''''(i)(4∆x)^4/4! + (e)f'''''(i)(4∆x)^5/5! + (e)f''''''(i)(4∆x)^6/6!+........

(g)f(i+5)= (g)f(i)+(g)f'(i)(5∆x)/1! + (g)f''(i)(5∆x)^2/2! + (g)f'''(i)(5∆x)^3/3! - (g)f''''(i)(5∆x)^4/4! + (g)f'''''(i)(5∆x)^5/5! + (g)f''''''(i)(5∆x)^6/6!+........

arranging taylor series expansion in taylor term:

	f(i)	∆xf(i)	∆x^2f''(i)	∆x^3f'''(i)	∆x^4f''''(i)	∆x^5f'''''(i)	∆x^6f''''''(i)
af(i)	a	0	0	0	0	0	0
(b)f(i+1)	b	b	b/2	b/6	b/24	b/120	b/720
(c)f(i+2)	c	2c	4c/2	8c/6	16c/24	32c/20	64c/720
(d)f(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120	729d/720
(e)f(i+4)	e	4e	16e/2	64e/6	256e/24	1024e/120	4096e/720
(g)f(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120	15625g/720
sum of above all	0	0	1	0	0	0	?

from above this we can write in matrix form AX=B

[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 1024/120 ; 0 1/120 32/120 243/120 1024/120 3125/120 ] [a b c d e g] = [ 0 0 1 0 0 0]

after solving the matrix the value we get:

a=3.75000

b=-12.83333

c= 17.83333

d=-13.00000

e=5.08333

g=-0.83333

dividing both the side by equation (1) by (∆x^2) we get

∂²u/∂x² = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5) %we can put the above vlues(a b c d e g) in it

putting the above all values we can get the equation.

USING SKEWED LEFT SIDED DIFFERENCE:

In this skewed left hand side difference method we have to define the five points on the left of (i).

the numerical stencil will be:

(i-5) (i-4) (i-3) (i-2) (i-1) (i)

the hypothesis eqquation will be

∂²u/∂x² = af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)

writing taylor series expansion for each term

(a)f(i-5)= (a)f(i)-(a)f'(i)(5∆x)/1! + (a)f''(i)(5∆x)^2/2! - (a)f'''(i)(5∆x)^3/3! + (a)f''''(i)(5∆x)^4/4! - (a)f'''''(i)(5∆x)^5/5! + (a)f''''''(i)(5∆x)^6/6!+.........

(b)f(i-4)=(b)f(i)-(b)f'(i)(4∆x)/1! + (b)f''(i)(4∆x)^2/2! - (b)f'''(i)(4∆x)^3/3! + (b)f''''(i)(4∆x)^4/4! - (b)f'''''(i)(4∆x)^5/5! + (b)f''''''(i)(4∆x)^6/6!+.......

(c)f(i-3)=(c)f(i)-(c)f'(i)(3∆x)/1! + (c)f''(i)(3∆x)^2/2! - (c)f'''(i)(3∆x)^3/3! + (c)f''''(i)(3∆x)^4/4! - (c)f'''''(i)(3∆x)^5/5! + (c)f''''''(i)(3∆x)^6/6!+..........

(d)f(i-2)= (d)f(i)-(d)f'(i)(2∆x)/1! + (d)f''(i)(2∆x)^2/2! - (d)f'''(i)(2∆x)^3/3! + (d)f''''(i)(2∆x)^4/4! - (d)f'''''(i)(2∆x)^5/5! + (d)f''''''(i)(2∆x)^6/6!+........

(e)f(i-1)= (e)f(i)-(e)f'(i)(∆x)/1! + (e)f''(i)(∆x)^2/2! - (e)f'''(i)(∆x)^3/3! + (e)f''''(i)(∆x)^4/4! - (e)f'''''(i)(∆x)^5/5! + (e)f''''''(i)(∆x)^6/6!+........

(g)f(i)=(g)f(i)

arranging taylor series expansion in taylor term:

	f(i)	∆xf(i)	∆x^2f''(i)	∆x^3f'''(i)	∆x^4f''''(i)	∆x^5f'''''(i)	∆x^6f''''''(i)
af(i-5)	a	-5a	25/a	-125a/6	625a/24	-3125a/120	15625a/720
(b)f(i-4)	b	-4b	8b	-64b/6	256b/24	-1024b/120	4096b/720
(c)f(i-3)	c	-3c	9c/2	-27c/6	81c/24	-243c/120	729c/120
(d)f(i-2)	d	-2d	2d	-8d/6	16d/24	32d/120	64d/720
(e)f(i-1)	e	-e	e/2	-e/6	e/24	-e/120	e/720
(g)f(i)	g	0	0	0	0	0	0
sum of above all	0	0	1	0	0	0	?

from above this we can write in matrix form AX=B

[1 1 1 1 1 1 ; -5 -4 -3 -2 -1 0 ; 25/2 8/2 9/2 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 ] [ a b c d e g]=[ 0 0 1 0 0 0]

solving the above matrix we get the values

a= -0.83333

b= 5.08333

c=-13.00000

d=17.83333

e=-12.83333

g=3.75000

dividing both the side by equation (1) by (∆x^2) we get

∂²u/∂x² = af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i) %we can put the above vlues(a b c d e g) in it

putting the above all values we can get the equation.

% function for central difference error

% creating the fuction for 4order approximation
function  difference_central_order =central_difference(x,dx)

%analytical function f(x)= exp(x)*cos(x);
% f''(x)=-2*exp(x)*sin(x);

analytical_deravative= -2*exp(x)*sin(x);
% values of coefficients from  taylor table
a=-0.083333;
b=1.333333;
c=-2.500000;
d=1.333333;
e=-0.83333;

%central differencing
%((a*f(x-2*dx))+(b*f(x-dx))+(c*f(x))+(d*f(x+dx))+(e*f(x+2*dx)))/dx^2;

central_differencing_fourth_order=((a*exp(x-2*dx)*cos(x-2*dx))+(b*exp(x-dx)*cos(x-dx))+(c*exp(x)*cos(x))+(d*exp(x+dx)*cos(x+dx))+(e*exp(x+2*dx)*cos(x+2*dx)))/dx^2;
difference_central_order= abs(central_differencing_fourth_order-analytical_deravative)
end

% function for skewed right side

% creating 4 order approximation using 2 order deravative by skewed right side difference
function error_skewed_right_side=skewed_right_side(x,dx)

%analytical function f(x)= exp(x)*cos(x);
% f''(x)=-2*exp(x)*sin(x);

analytical_deravative= -2*exp(x)*sin(x);
% values of coefficients from  taylor table
a=3.75000;
b=-12.83333;
c= 17.83333;
d=-13.00000;
e=5.08333;
g=-0.83333;

skewed_right_side_fourth_order=((a*exp(x)*cos(x))+(b*exp(x+dx)*cos(x+dx))+(c*exp(x+2*dx)*cos(x+2*dx))+(d*exp(x+3*dx)*cos(x+3*dx))+(e*exp(x+4*dx)*cos(x+4*dx))+(g*exp(x+5*dx)*cos(x+5*dx)))/dx^2;

error_skewed_right_side=abs(skewed_right_side_fourth_order-analytical_deravative)
end

%function for skewed left side

 % creating 4 order approximation using 2 order deravative by skewed left side difference
function error_skewed_left_sided= skewed_left_side(x,dx)
%analytical function f(x)= exp(x)*cos(x);
% f''(x)=-2*exp(x)*sin(x);

analytical_deravative= -2*exp(x)*sin(x);
% values of coefficients from  taylor table
a= -0.83333;
b= 5.08333;
c=-13.00000;
d=17.83333;
e=-12.83333;
g=3.75000;

skewed_left_side_fourth_order=((a*exp(x-5*dx)*cos(x-5*dx))+(b*exp(x-4*dx)*cos(x-4*dx))+(c*exp(x-3*dx)*cos(x-3*dx))+(d*exp(x-2*dx)*cos(x-2*dx))+(e*exp(x-dx)*cos(x-dx))+(g*exp(x)*cos(x)))/dx^2;
error_skewed_left_sided=abs(skewed_left_side_fourth_order-analytical_deravative);
end

%programing

clear all
clc
close all
 x=pi/3;
 dx=linspace(pi/4,pi/4000,30);

 for i=1:length(dx)
    central_difference_error (i) =central_difference(x,dx(i));
     skewed_right_side_error (i) =skewed_right_side(x,dx(i));
     skewed_left_side_error (i) =skewed_left_side(x,dx(i));
 end
 % plotting graph of dx vs error
 figure(1)
 plot(dx,central_difference_error,'LineWidth',2)
 hold on
 plot(dx,skewed_right_side_error,"Color",'r','LineWidth',2)
 hold on
 grid on
 plot(dx,skewed_left_side_error,"Color",'g','LineWidth',2)
 xlabel('dx')
 ylabel('error')
 legend('central difference error','skewed right sided','skewed left side')

 figure(2)
  loglog(dx,central_difference_error,'LineWidth',2)
 hold on
 loglog(dx,skewed_right_side_error,"Color",'r','LineWidth',2)
 hold on
 grid on
 loglog(dx,skewed_left_side_error,"Color",'g','LineWidth',2)
 xlabel('dx')
 ylabel('error')
 legend('central difference error','skewed right sided','skewed left side')

%result graph

here we can easily notice central difference approximation produce more accurate solution.