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

1)Central Difference Fourth Order Approximation:-

`(del^2u)/(del^2x)cong a[f(i-2)]+b[f(i-1)]+c[f(i)]+d[f(i+1)]+e[f(i+2)]`

	f(i)	`Delta(x)f'(x)`	`Delta(x)f''(x)`	`Delta(x)f'''(x)`	`Delta(x)f^(iv)(x)`
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`=`(del^2u)/(delx^2)`	0	0	1	?	?

The solutions are:-

2)Right Skewed Fourth Order Approximation:-

`(del^2u)/(del^2x)cong a[f(i)]+b[f(i+1)]+c[f(i+2)]+d[f(i+3)]+e[f(i+4)+g[f(i+5)]`

	f(i)	`Delta(x)f'(x)`	`Delta(x)f''(x)`	`Delta(x)f'''(x)`	`Delta(x) f^(iv)(x)`	`Delta(x)f^(v)(x)`
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	8c/6	16c/24	32c/120
df(i+3)	d	3d	9d/2	27d/6	81d/24	243d/120
ef(i+4)	e	4e	8e	64e/6	256e/24	1024e/120
gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120
`sum`=`(del^2u) /(delx^2)`	0	0	1	0	0	?

The solutions are:-

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

3)Left Skewed Fourth Order Approximation:-

`(del^2u)/(del^2x)cong a[f(i-5)]+b[f(i-4)]+c[f(i-3)]+d[f(i-2)]+e[f(i-1)+g[f(i)]`

	f(i)	`Delta(x)f'(x)`	`Delta(x)f''(x)`	`Delta(x)f'''(x)`	`Delta(x) f^(iv)(x)`	`Delta(x)f^(v)(x)`
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	-8c/6	16c/24	-32c/120
df(i+3)	d	-3d	9d/2	-27d/6	81d/24	-243d/120
ef(i+4)	e	-4e	8e	-64e/6	256e/24	-1024e/120
gf(i+5)	g	-5g	25g/2	-125g/6	625g/24	-3125g/120
`sum`=`(del^2u) /(delx^2)`	0	0	1	0	0	?

The solutions are:-

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

clear all
close all
clc

%Central Difference (i-2,i-1,1,i+1,i+2)
   C = [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];
   A = inv(C)*B

% Skewed Right Side Difference (i,i+1,1,i+2,i+3,i+4,i+5)
  k = [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];
  l = [0; 0; 1; 0; 0; 0];
  m = inv(k)*l

% Skewed Left Side Difference (i,i-1,i-2,i-3,i-4,i-5)
  u = [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];
  v = [0; 0; 1; 0; 0; 0];
  w = inv(u)*v
% Taylor Table Formulation: 

% Function - f(x) = e^x(cos(x))                                           
%            f'(x) = e^x(cos(x) - sin(x))                                 
%            f''(x) = -2e^x*sin(x)                                        
close all
clc

% Input Variables
x = pi/3;
dx = linspace(pi/40,pi/400,100); % step size

% Iterative Loop 
for i = 1: length(dx)
     [central_diff(i),skewed_right(i),skewed_left(i),central_error(i),skewed_right_error(i),skewed_left_error(i)] = taylor_scheme(x,dx(i))
end


% Plot                                          


figure(1)
loglog(dx,skewed_left_error,'r')
hold on
loglog(dx,skewed_right_error,'b')
hold on
loglog(dx,central_error,'g')
legend('Skewed Right-Sided Difference','Skewed Left-Sided Difference','Central Difference')
xlabel('Length of Increment dx')
ylabel('Truncation Error')


%FUNCTIONS                          


function [central_diff,skewed_right_diff,skewed_left_diff,central_error,skewed_right_error,skewed_left_error] = taylor_scheme(x,dx)

% f(x) = exp(x)*cos(x)
% simplify the equation 
  k1 = exp(x)*cos(x);                   
  k2 = exp(x+dx)*cos(x+dx);              
  k3 = exp(x+(2*dx))*cos(x+(2*dx));      
  k4 = exp(x+(3*dx))*cos(x+(3*dx));      
  k5 = exp(x+(4*dx))*cos(x+(4*dx));      
  k6 = exp(x-(4*dx))*cos(x-(4*dx));     
  k7 = exp(x-(3*dx))*cos(x-(3*dx));      
  k8 = exp(x-(2*dx))*cos(x-(2*dx));      
  k9 = exp(x-dx)*cos(x-dx);              
  k10 = exp(x-5*dx)*cos(x-5*dx);         
  k11 = exp(x+5*dx)*cos(x+5*dx);        

% Analytical Derivative
    exact_derivative = -2*exp(x)*sin(x);

% Central Difference 
% CD = (af(i-2) + bf(i-1) + cf(i) + df(i+1) + ef(i+2))/(dx^2)
   c = [-1/12; 4/3; -5/2; 4/3; -1/12];
   central_diff = (c(1)*k8 + c(2)*k9 + c(3)*k1 + c(4)*k2 + c(5)*k3)/(dx^2);

% Skewed Right Side Difference 
% SRD = (af(i) + bf(i+1) + cf(i+2) + df(i+3) + ef(i+4) + gf(i+5))/(dx^2)
  m = [15/4; -77/6; 107/6; -26/2; 61/12; - 5/6];
  skewed_right_diff = (m(1)*k1 + m(2)*k2 + m(3)*k3 + m(4)*k4 + m(5)*k5 + m(6)*k11)/(dx^2);

% Skewed Left Side Difference 
  % SLD = (af(i) + bf(i-1) + cf(i-2) + df(i-3) + ef(i-4) + gf(i-5))/(dx^2)
  w = [15/4; -77/6; 107/6; -26/2; 61/12; - 5/6];
  skewed_left_diff = (w(1)*k1 + w(2)*k9 + w(3)*k8 +w(4)*k7 +w(5)*k6 +w(6)*k10)/(dx^2);

 % Error 
  central_error = abs(central_diff - exact_derivative);
  skewed_right_error = abs(skewed_right_diff - exact_derivative);
  skewed_left_error = abs(skewed_left_diff - exact_derivative);

end

Plot:-

Conclusion:-

From the above result, we can observe that the Truncation error for different methods is as below

Central Difference < Skewed_Right_Sided_Difference < Skewed_Left_Sided_Difference.