clear all
close all
clc
x = 2;
% To determine the Analytical solution of f(x) = exp(x)*cos(x)
% that is 1st order derivative as exp(x)*(-sin(x)) and
%2nd order derivative as exp(x)*(-cos(x)) and
ana = 2*exp(x)*(-sin(x));
%Initial inputs from hand calculations and taylor table
dx = 0.2; %Element size
%Let[Ac] be matrix containing a,b,c,d and e for
%centralised differentiation of 2nd order Differential eqn.
%Let[Ar] be matrix containing a,b,c,d,e and f for
%Right skewed differentiation of 2nd order Differential eqn.
%Let[Al] be matrix containing a,b,c,d,e and f for
%Left skewed differentiation of 2nd order Differential eqn.
%Let [C1] be matrix containing co-efficients of central differencing
%Let [C2] be matrix containing co-efficients of right skewed differencing
%Let [C3] be matrix containing co-efficients of left skewed differencing
%Let [R1] be matrix RHS of central differencing
%Let [R2] be matrix RHS of right skewed differencing
%Let [R3] be matrix RHS of left skewed differencing
C1 = [1 1 1 1 1
-2*(dx) -dx 0 dx 2*(dx)
2*(dx^2) 0.5*(dx^2) 0 0.5*(dx^2) 2*(dx^2)
-(4/3)*(dx^3) (dx^3)/6 0 (dx^3)/6 (4/3)*(dx^3)
(16/24)*(dx^4) (dx^4)/24 0 (dx^4)/24 (16/24)*(dx^4)];
R1 = [0
0
1
0
0];
%Define arbitrary co-efficient Y1 to determine inverse of C1
Y1= inv(C1);
A1 = Y1 * R1;
% Also
C2 = [1 1 1 1 1 1
0 dx 2*dx 3*dx 4*dx 5*dx
0 (dx^2)/2 ((2*dx)^2)/2 ((3*dx)^2)/2 ((4*dx)^2)/2 ((5*dx)^2)/2
0 (dx^3)/6 ((2*dx)^3)/6 ((3*dx)^3)/6 ((4*dx)^3)/6 ((5*dx)^3)/6
0 (dx^4)/24 ((2*dx)^4)/24 ((3*dx)^4)/24 ((4*dx)^4)/24 ((5*dx)^4)/24
0 (dx^5)/120 ((2*dx)^5)/120 ((3*dx)^5)/120 ((4*dx)^5)/120 ((5*dx)^5)/120];
R2 = [0
0
1
0
0
0];
%Define arbitrary co-efficient Y2 to determine inverse of C2
Y2= inv(C2);
Ar = Y2 * R2;
% Also
C3 = [1 1 1 1 1 1
-5*dx -4*dx -3*dx -2*dx -dx 0
((5*dx)^2)/2 ((4*dx)^2)/2 ((3*dx)^2)/2 ((2*dx)^2)/2 (dx^2)/2 0
-((5*dx)^3)/6 -((4*dx)^3)/6 -((3*dx)^3)/6 -((2*dx)^3)/6 -(dx^3)/6 0
((5*dx)^4)/24 ((4*dx)^4)/24 ((3*dx)^4)/24 ((2*dx)^4)/24 (dx^4)/24 0
-((5*dx)^5)/120 -((4*dx)^5)/120 -((3*dx)^5)/120 -((2*dx)^5)/120 -(dx^5)/120 0];
R3 = [0
0
1
0
0
0];
%Define arbitrary co-efficient Y2 to determine inverse of C2
Y3= inv(C3);
Al = Y3 * R3;
%To calculate 2nd order derivative by Taylor table method
% that is,
%By Center Difference
%f"(x) = a*f(x-2*dx) + b *f(x-dx) + c*f(x) + d*f(x+dx) + e*f(x+2dx)
%Define fc to determine numerical value of f(x) by center differencing
%Also define f(x-2*dx) f(x-dx) f(x) f(x+dx)f(x+2dx)
%to determine numerical value of f(x) by center differencing
%in matrix fdc
fdc = [exp(x-(2*dx))*cos(x-(2*dx)) exp(x-dx)*cos(x-dx) exp(x)*cos(x) exp(x+dx)*cos(x+dx) exp(x+(2*dx))*cos(x+(2*dx))];
fc = fdc * A1;
%Define error_c to calculate the error
error_c = abs(fc-ana);
%By Right skewed Difference
%f"(x) = a*f(x) + b *f(x+dx) + c*f(x+2*dx) + d*f(x+3dx) + e*f(x+4dx) + f*f(x+5dx)
%Define fr to determine numerical value of f(x) by right skewed differencing
%Also define fc1 =f(x) fc2= f(x+dx) fc3= f(x+3dx) fc4 = f(x+4*dx) fc5 =f(x+5dx)
%in matrix fdr
fdr = [exp(x)*cos(x) exp(x+dx)*cos(x+dx) exp(x+(2*dx))*cos(x+(2*dx)) exp(x+(3*dx))*cos(x+(3*dx)) exp(x+(4*dx))*cos(x+(4*dx)) exp(x+(5*dx))*cos(x+(5*dx))];
fr = fdr * Ar;
%Define error_r to calculate the error
error_r = abs(fr-ana);
%By Left skewed Difference
%f"(x) = a*f(x-4*dx) + b *f(x-3*dx) + c*f(x-2*dx) + d*f(x-dx) + e*f(x)
%Define fl to determine numerical value of f(x) by left skewed differencing
%By Right skewed Difference
%f"(x) = a*f(x) + b *f(x+dx) + c*f(x+2*dx) + d*f(x+3dx) + e*f(x+4dx) + f*f(x+5dx)
%Define fr to determine numerical value of f(x) by right skewed differencing
%Also define f(x-5dx) f(x-4dx) f(x-3dx) f(x-2*dx) f(x-dx) f(x)
%in matrix fdl
fdl = [exp(x-(5*dx))*cos(x-(5*dx)) exp(x-(4*dx))*cos(x-(4*dx)) exp(x-(3*dx))*cos(x-(3*dx)) exp(x-(2*dx))*cos(x-(2*dx)) exp(x-dx)*cos(x-dx) exp(x)*cos(x)];
fl = fdl * Al;
%Define error_l to calculate the error
error_l = abs(fl-ana);
e = [error_c error_r error_l];
%Plotting a bar chart by titling every bar corresponding
%to their respective errors
%Blue for error in 1 order approximation
%Yellow for error in 2nd order appproximation
%Green for error in 4th order approximation
x = categorical({'Center differenc error' , 'Skewed right difference error','Skewed left difference error'});
x = reordercats(x,{'Center differenc error' , 'Skewed right difference error','Skewed left difference error'});
y = [ error_c error_r error_l] ;
figure
b = bar (x,y)
%To define different color codes for different errors
b.FaceColor = 'flat'
b.CData(2,:) = [1 8 .3]
b.CData (3,:) = [0.2 0.6 .5]
% To define the title of chart and axes
title ('Error comparision')
xlabel ('Order of differencing')
ylabel ('Error')
figure
plot(x,y)
title ('Error v/s Order of approximation')
xlabel('Order of differencing')
ylabel('Error')