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

%%---- 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')

CODE EXPLANATION:

• Each scheme is coded in the form of 'function'.
• The equation for all the schemes were solved initially on paper and in order to obtain the coefficients of the equation, they were solved using MATLAB (you can see in code).
• It will ask user to input the value of x and range for dx respectively and these values are supplied to the function. After solving the function, the function itself is called in the main code.
• The main code will give the comparison of the errors of all the schemes.

THEORY EXPLANATION:

• The idea behinde using central, right hand skewed and left hand skewed schemes is to obtiain higher order approximation at the nodes which do not have left or right point (they are starting point).
• Therefore, one can get the higher order approximation using right hand and left hand scheme.

HOW TO DECIDE NODES IN TAYLOR SERIES:

• For skewed scheme, Number of nodes = order of derivative + order of approximation. In our case, order of derivative is 2nd order and order of approximation is 4th order. Therefore, Number of node = 2+4=6.
• For central scheme, Number of nodes = order of derivative + order of approximation -1.In our case, order of derivative is 2nd order and order of approximation is 4th order. Therefore, Number of node = 2+4-1=5.
• For right hand skew, [af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)]/(dx^2) (in our case)
• For left hand skew, [af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)]/(dx^2) (in our case)
• For central scheme, [af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)]/(dx^2) (in our case)
• Initially, one must calculate the coefficients and then solve the equations.