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

% Main function
clc
clear all
close all


x = pi/30;
dx= linspace(pi/40, pi/4000,20);

for i = 1:length(dx)

     error_cd(i) = central_diff (x,dx(i))

     error_srd(i) = right_diff (x,dx(i))

     error_sld(i) = left_diff (x,dx(i))


end

 loglog(dx, error_cd, 'r', dx, error_srd, 'g', dx, error_sld, 'y')
 title ('Fourth order approximation')

 grind on
 
 x label ('dx');
 y label ('Errors dx');
 legend ('error_cd', 'error_srd', 'error_sld')

% Right diff

%Right skewed scheme 

function error_srd = right_diff (x,dx)

%Analytical function f(x) = exp(x)*cos(x);
%SECOND ORDER DERIVATIVE f"(x) = -2*exp(x).sin(x);

analytical_derivative = -2*exp(x)*sin(x);

%Coefficient values

a= 3.750000;
b= -12.833333;
c= 17.833333;
d= -13.000000;
e= 5.083333;
g = -0.833333;

%right skewed differencing expansion -

% ((a*f(x)+(b*f(x+dx))+(c*f(x+2dx))+(d*f(x+3dx)+(e*f(x+4*dx)))/(dx^2);

rightskewed_diff_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_srd = abs (rightskewed_diff_fourth_order - analytical_derivative);

end

% Left diff


% Skewed left scheme

function error_sld = left_diff (x,dx)


%Analytical function f(x) = exp(x)*cos(x);
%SECOND ORDER DERIVATIVE f"(x) = -2*exp(x).sin(x);

analytical_derivative = -2*exp(x)*sin(x);

%Coefficient values

a= -0.833333;
b= 5.083333;
c= -13.000000;
d= 17.833333;
e= -12.833333;
g = 3.750000;

%left skewed differencing expansion -

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

leftskewed_diff_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_sld = abs (leftskewed_diff_fourth_order - analytical_derivative);

end


% Central diff
% CENTRAL DIFFERENCE

function error_cd = central_diff (x,dx)


%Analytical function f(x) = exp(x)*cos(x);
%SECOND ORDER DERIVATIVE f"(x) = -2*exp(x).sin(x);

analytical_derivative = -2*exp(x)*sin(x);

%Coefficient values

a=-0.083333;
b=1.333333;
c=-2.500000;
d=1.333333;
e=-0.083333;

%central differencing expansion -

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

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

error_cd = abs (central_diff_fourth_order - analytical_derivative);


end

Derivation of 4th order approximations of the second order derivative for -

1. Central difference

Consider (i-2), (i-1), (i), (i+1) and (i+2) as the nodes in the numerical stencil.

`(del^2u)/(delx^2) = a.f(i-2)+b.f(i-1)+c.f(i)+d.f(i+1)+e.f(i+2)` _____ (1)

Upon expansion of terms and arrangement, we get

	f(i)	`∆`x.f'(i)`	`∆^2`x.f''(i)`	`∆^3`x.f'''(i)`	`∆^4`x.f''''(i)`	`∆^5`x.f'''''(i)`	`∆^6`x.f''''''(i)`
a.f(i-2)	a	-2a	4a/2	-8a/6	16a/24	-32a/120	64a/120
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
`(del^2u)/(delx^2)`	0	0	1	0	0	-	-

The equations thus formed are as follows.

a+b+c+d+e=0

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

2a+(b/2)+0+(d/2)+2e = 1

-(8a/6)-(b/6)+0+(d/6)+(8e/6) = 0

(16a/24)+(b/24)+0+(d/24)+(16e/24) = 0

Then matrix 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]

On doing matrix multiplication and solving - a = -0.0833, b = 1.3333, c = -2.5000, d = 1.3333, e = -0.0833

Here a = e and b = d.

Similarly, for

2. Right skewed scheme -

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

3. Left skewed scheme -

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