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

Aim:

To derive the following 4th order approximations of the second-order derivative.

1. Central difference

2. Skewed right-sided difference

3. Skewed left-sided difference

%%---- CENTRAL SCHEME ----%%

% Function for 4th ORDER APPROX OF 2nd ORDER DERIVATIVE 
% function = exp(x)*cos(x)
% analytical_derivative = - exp(x)*sin(x) + exp(x)*cos(x)

% CODE
function output = four_two_central(x,dx) % central 
% For central scheme coefficient calculation
a1 =[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];
b1=[0; 0; 1; 0; 0];
coeff= linsolve(a1,b1);

% Analytical Solution
% analytical_derivative = (exp(x) * (cos(x)-sin(x)));
 analytical_derivative = -2*(exp(x) * (sin(x)));
 
% Numerical Solution
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))));
fourth_numerical_central = (((coeff(1))*a)+((coeff(2))*b)+((coeff(3))*c)+((coeff(4))*d)+(((coeff(5))*e)))/(dx.^2);

% Error 
output = abs(fourth_numerical_central - analytical_derivative);
end

%%---- RIGHT HAND SCHEME ----%%

% Function for 4th ORDER APPROX OF 2nd ORDER DERIVATIVE 
% function = exp(x)*cos(x)
% analytical_derivative = - exp(x)*sin(x) + exp(x)*cos(x)

% CODE
function output = fourth_forward_skew(x,dx) % (right hand side)
% For forward skewed scheme coefficient calculation (right hand side)
c1 = [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];
d1 = [0; 0; 1; 0; 0; 0];
coeff= linsolve(c1,d1);

% Analytical Solution
% analytical_derivative = (exp(x) * (cos(x)-sin(x)));
 analytical_derivative = -2*(exp(x) * (sin(x)));
 
 % Numerical Solution
p= (exp(x)*cos(x));
 q= ((exp(x+dx)).*(cos(x+dx)));
 r= ((exp(x+(2*dx))).*(cos(x+(2*dx))));
 s= ((exp(x+(3*dx))).*(cos(x+(3*dx))));
 t= ((exp(x+(4*dx))).*(cos(x+(4*dx))));
 u= ((exp(x+(5*dx))).*(cos(x+(5*dx)))); 
forward_skew = (((15/4)*p)+((-77/6)*(q))+((107/6)*(r))+((-13)*(s))+((61/12)*(t))+((-5/6)*(u)))/(dx.^2);

% Error 
output = abs(forward_skew - analytical_derivative);
end

%%---- LEFT HAND SCHEME---- %%

% Function for 4th ORDER APPROX OF 2nd ORDER DERIVATIVE 
% function = exp(x)*cos(x)
% analytical_derivative = - exp(x)*sin(x) + exp(x)*cos(x)

% CODE
function output = fourth_backward_skew(x,dx) % Left hand side
% For backward skewed scheme coefficient calculation (left hand side)
l1 = [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];
m1 = [0; 0; 1; 0; 0; 0];
coeff= linsolve(l1,m1);
% Analytical Solution
% analytical_derivative = (exp(x) * (cos(x)-sin(x)));
  analytical_derivative = -2*(exp(x) * (sin(x)));
 
 % Numerical Solution
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)));
backward_skew = ((-5/6*(a))+(61/12*(b))+(-13*(c))+(107/6*(d))+(-77/6*(e))+(15/4*(g)))/(dx.^2);

% Error
output = abs(backward_skew - analytical_derivative);
end

%%---- MAIN CODE ----%%

% Main Code
n = input('Specify Value of x: ');
m = input('Input start value of dx: ');
o = input('Input end value of dx: ');
p = input('Input number of points: ');
x= n;
dx= linspace(m,o,p);

% Analytical Solution
analytical_derivative = (exp(x) * (cos(x)-sin(x)));

% For loop for ERROR calculation
for i = 1:length(dx)
   central_error(i) = four_two_central(x,dx(i)); 
   forward_skew_error(i) = fourth_forward_skew(x,dx(i));
   backward_skew_error(i) = fourth_backward_skew(x,dx(i));
end

% Plotting
figure(1)
loglog(dx,central_error,'r')
hold on
loglog(dx,forward_skew_error,'g')
hold on
loglog(dx,backward_skew_error,'b')
xlabel('dx')
ylabel('Error')
grid on;
legend('Central Error','Forward Error','Backward Error')