Week 3 - Taylor table method and Matlab code

function Y = skewed_right_difference(C,D)

Y=C\D;
end

function Z = skewed_left_difference(E,F)

Z=E\F;
end

function X = Central_difference(A,B)

X=A\B;
end

% Function for central differencing scheme for 2nd order
function out=central_difference_approximation(x,dx,X)
  % Analytical function f(x)=exp(x)*cos(x)
  %Analytical derivative
  %f''(x)=-2*exp(x)*sin(x)
  analytical_derivative =(-2*exp(x)*sin(x));
  % Numerical derivative
  %Fourth order approximation

 % Central differencing=[a*f(x-2dx)+b*f(x-dx)+c*f(x)+d*f(x+dx)+e*f(x+2dx)]/dx^2
 Central_Differencing=(X(1)*exp(x-(2*dx))*cos(x-(2*dx))+X(2)*exp(x-(1*dx))*cos(x-(1*dx))+X(3)*exp(x)*cos(x)+X(4)*exp(x+(1*dx))*cos(x+(1*dx))+X(5)*exp(x+(2*dx))*cos(x+(2*dx)))/(dx^2);
 % Error
 out=abs(Central_Differencing-analytical_derivative);
end

%Function for skewed right side differencing scheme
function out=skewed_rightside_approximation(x,dx,Y)
   % f(x)=exp(x)*cos(x)
   analytical_derivative=-2*exp(x)*sin(x);
   %Numerical Derivative
   %Fourth order approximation
  
   % skewed right side differencing=[a*f(x)+b*f(x+dx)+c*f(x+2dx)+d*f(x+3dx)+e*f(x+4dx)+g*f(x+5dx)]/(dx^2)
   skewed_rightsided_differencing=(Y(1)*exp(x)*cos(x))+(Y(2)*exp(x+dx)*cos(x+dx))+(Y(3)*exp(x+2*dx)*cos(x+2*dx))+(Y(4)*exp(x+3*dx)*cos(x+3*dx))+(Y(5)*exp(x+4*dx)*cos(x+4*dx))+(Y(6)*exp(x+5*dx)*cos(x+5*dx))/(dx^2);
   out=abs(skewed_rightsided_differencing-analytical_derivative);
end
   

%Function for skewed left side differencing scheme
function out=skewed_leftside_approximation(x,dx,Z)
   % f(x)=exp(x)*cos(x)
   analytical_derivative=-2*exp(x)*sin(x);
   %Numerical Derivative
   %Fourth order approximation
  
   % skewed left side differencing=[a*f(x)+b*f(x-dx)+c*f(x-2dx)+d*f(x-3dx)+e*f(x-4dx)+g*f(x-5dx)]/(dx^2)
   skewed_leftsided_differencing=(Z(1)*exp(x)*cos(x))+(Z(2)*exp(x-dx)*cos(x-dx))+(Z(3)*exp(x-2*dx)*cos(x-2*dx))+(Z(4)*exp(x-3*dx)*cos(x-3*dx))+(Z(5)*exp(x-4*dx)*cos(x-4*dx))+(Z(6)*exp(x-5*dx)*cos(x-5*dx))/(dx^2);
   out=abs(skewed_leftsided_differencing-analytical_derivative);
end

%Numerical solution 
clear all
close all
clc
% Analytical function to be assumed as exp(x)*cos(x)
%f''(x)=-2exp(x)*sin(x)

%Assuming suitable value of x and dx
x=pi/3;
dx=linspace(pi/4,pi/4000,100);
central_error = zeros(1,length(dx));
skewed_right_error=zeros(1,length(dx));
skewed_left_error=zeros(1,length(dx));
A=[1 1 1 1 1;-2 -1 0 1 2;4/2 1/2 0 1/2 4/2;-8/6 -1/6 0 1/6 8/6;16/24 1/24 0 1/24 16/24];
B=[0; 0; 1; 0; 0];
X = Central_difference(A,B);
C=[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];
D=[0; 0; 1; 0; 0;0];
Y = skewed_right_difference(C,D);

E=[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];
F=[0; 0; 1; 0; 0;0];
Z = skewed_left_difference(E,F);
%loop to compute the numerical solution for each value of dx
for i=1:length(dx)
    central_error(i)=central_difference_approximation(x,dx(i),X);
    skewed_right_error(i)=skewed_rightside_approximation(x,dx(i),Y);
    skewed_left_error(i)=skewed_leftside_approximation(x,dx(i),Z);
     %FDS_error=forward_approximation(x,dx(i));
     %BDS_error=backward_approximation(x,dx(i));
end
% Plotting the graphs and comparing the error
figure(1)
loglog(dx,central_error,'r','linewidth',2)
hold on
loglog(dx,skewed_right_error,'g','linewidth',2)
loglog(dx,skewed_left_error,'b','linewidth',2)
xlabel('grid size(dx)')
ylabel('Absolute Error')
legend('central error','skewed right error','skewed left error')
title('variation of different error')