Taylor table using matlab code
CENTRAL DIFFERENCING SCHEME:
Since CDS uses points on either side, the equation is as follows
`(del^2f)/(delx^2)~= a.f(i-2)+ b.f(i-1)+c.f(i)+d.f(i+1)+e.f(i+2)`
(assuming this is equation 1)
The taylor table is tabulated below
f(i) f'(i) f''(i) f'''(i) f''''(i) f'''''(i) f''''''(i)
a.f(i-2) a -2a 4a/2! -8a/3! 16a/4! 32a/5! 64a/6!
b.f(i-1) b -b b -b/3! b/4! -b/5! b/6!
c.f(i) c 0 0 0 0 0 0
d.f(i+1) d d d/2! d/3! d/4! d/5! d/6!
e.f(i+2) e 2e 4e/2! 8e/3! 16e/4! 32e/5! 64e/6!
Summing them up and equating the coefficients of f(i), f'(i),f'''(i) and f''''(i) to 0 and f''(i) to 1,we get a system of linear equations in the following matrix form
`[[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]][[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`
Using matlab we get
a =-1/12; b=4/3 ; c=-5/2; d=4/3 ; e -1/12
`(del^2f)/(delx^2)~=(-1/12.f(i-2)+4/3f(i-1)-5/2f(i)+4/3f(i+1)-1/12f(i+2))/(dx)^2 +O(dx)^4`
SKEWED RIGHT SIDED DIFFERENCE SCHEME (FDS)
The approximation is given by
`(del^2f)/(delx^2)~=a.f(i)+b.f(i+1)+c.f(i+2)+d.f(i+3)+e.f(i+4)+h.f(i+5)`
(Assuming this is equation 2)
The taylor table is given by
f(i) f'(i) f''(i) f'''(i) f''''(i) f'''''(i) f''''''(i)
a.f(i) a 0 0 0 0 0 0
b.f(i+1) b b b/2! b/3! b/4! b/5! b/6!
c.f(i+2) c 2c 4c/2 8c/3! 16c/4! 32c/5! 64c/6!
d.f(i+3) d 3d 9d/2! 27d/3! 81d/4! 243d/5! 729d/6!
e.f(i+4) e 4e 16e/2! 32e/3! 64e/4! 256e/5! 1024e/6!
h.f(i+5) h 5h 25h/2! 125h/3! 625h/4! 3125h/5! 15625h/6!
Summing it up and equating the coefficients of f(i), f'(i),f'''(i) ,f''''(i) and f'''''(i) to 0 and f''(i) to 1, we get the following system of linear equations in matrix form
`[[1,1,1,1,1],[0,1,2,3,4],[0,1/2,2,9/2,8],[0,1/6,4/3,9/2,32/3],[0,1/24,2/3,27/8,32/3],[0,1/120,32/120,243/120,1024/120]][[a],[b],[c],[d],[e],[h]]=[[0],[0],[1],[0],[0],[0]]`
Using matlab , we get
a=15/4; b=-77/6;c= 107/6; d=-13; e= 61/12;h= -5/6
substituting these values in equation 2, we get
`(del^2f)/(delf^2)~=(15/4f(i)-77/6f(i+1)+107/6f(i+2)-13f(i+3)+61/12f(i+4)-5/6f(i+5))/(dx)^2+O(dx)^4`
SKEWED LEFT SIDED DIFFERENCE SCHEME:
The approximation is given as follows
`(del^2f)/(delx^2)~= a.f(i-5)+b.f(i-4)+c.f(i-3)+d.f(i-2)+e.f(i-1)+h.f(i)`
(assuming this is equation 3)
Using taylor table method
f(i) f'(i) f''(i) f'''(i) f''''(i) f'''''(i) f''''''(i)
a.f(i-5) a -5a 25a/2! -125a/3! 625a/4! -3125a/5! 15625a/6!
b.f(i-4) b -4b 16b/2! -64b/3! 256b/4! -1024b/5! 4096/6!
c.f(i-3) c -3c 9c/2 ! -27c/3! 81c/4! -243c/5! 729c/6!
d.f(i-2) d -2d 4d/2! -8d/3! 16d/4! 32d/5! 64d/6!
e.f(i-1) e -e e/2! -e/3! e/4! -e/5! e/6!
h.f(i) 0 0 0 0 0 0 0
Summing it up and equating the coefficients of f(i),f'(i),f'''(i),f''''(i) ,f'''''(i) to 0 and f''(i) to 1 ,we get the following system of equations in matrix form
`[[1,1,1,1,1,1],[-5,-4,-3,-2,-1,0],[25/2,16/2,9/2,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]] [[a],[b],[c],[d],[e],[h]]=[[0],[0],[1],[0],[0],[0]]`
Using Matlab , we get
a=-5/6; b=61/12; c=-13; d=107/6; e=-77/6; h=15/4
substituting these values in equation 3,
`(del^2f)/(delx^2)~=(-5/6f(i-5)+61/12f(i-4)-13f(i-3)+107/6f(i-2)-77/6f(i-1)+15/4f(i))/dx^2+O(dx)^4`
MAIN CODE
clear all
close all
clc
% analytical function given = exp(x)*cos(x)
%analytical derivative= -2exp(x)*sin(x)
x= pi/2;
dx= linspace(0.01,0.1,10);
for i= 1:length(dx)
analytical_derivative= -2*exp(x)*(sin(x));
CDS_fourth_order_approx(i)= CDS(x,dx(i));
FDS_fourth_order_approx(i)= FDS(x,dx(i));
BDS_fourth_order_approx(i)= BDS(x,dx(i));
error_CDS(i)=abs(CDS_fourth_order_approx(i)-analytical_derivative);
error_FDS(i)=abs(FDS_fourth_order_approx(i)-analytical_derivative);
error_BDS(i)=abs(BDS_fourth_order_approx(i)-analytical_derivative);
end
%plots
%Figure(1)
loglog(dx,error_CDS,'r','marker','*','linewidth',1);
hold on
loglog(dx,error_FDS,'b','marker','square','linewidth',1);
loglog(dx,error_BDS,'g','marker','diamond','linewidth',1);
grid on
xlabel('grid size')
ylabel('error')
legend('central differencing scheme','forward differencing scheme','backward differencing scheme');
saveas(gca,'taylor_table_comparison.jpg')
% slope
log_error_CDS=log10(error_CDS);
log_error_FDS=log10(error_FDS);
log_error_BDS=log10(error_BDS);
log_dx=log10(dx);
slope_error_CDS=gradient(log_error_CDS(5),log_error_CDS(6),log_dx(5),log_dx(6));
slope_error_FDS=gradient(log_error_FDS(5),log_error_FDS(6),log_dx(5),log_dx(6));
slope_error_BDS=gradient(log_error_BDS(5),log_error_FDS(6),log_dx(5),log_dx(6));
% printing the slope
fprintf('slope using CDS is :%g\n',slope_error_CDS);
fprintf('slope using FDS is :%g\n',slope_error_FDS);
fprintf('slope using BDS is :%g\n',slope_error_BDS);
% function for CDS
function CDS_fourth_order = CDS(x,dx)
a=-1/12;
b=4/3;
c=-5/2;
d=4/3;
e=-1/12;
f_i_minus_2= e_(x-2*dx);
f_i_minus_1= e_(x-dx);
f_i= e_(x);
f_i_plus_1= e_(x+dx);
f_i_plus_2= e_(x+2*dx);
CDS_fourth_order=((a*f_i_minus_2)+(b*f_i_minus_1)+(c*f_i)+(d*f_i_plus_1)+(e*f_i_plus_2))/(dx^2);
end
% function for FDS
function FDS_fourth_order= FDS(x,dx)
a=15/4;
b=-77/6;
c=107/6;
d=-13;
e=61/12;
h=-5/6;
f_i=e_(x);
f_i_plus_1=e_(x+dx);
f_i_plus_2=e_(x+2*dx);
f_i_plus_3=e_(x+3*dx);
f_i_plus_4=e_(x+4*dx);
f_i_plus_5=e_(x+5*dx);
FDS_fourth_order= ((a*f_i)+(b*f_i_plus_1)+(c*f_i_plus_2)+(d*f_i_plus_3)+(e*f_i_plus_4)+(h*f_i_plus_5))/(dx^2);
end
% function for BDS
function BDS_fourth_order= BDS(x,dx)
a=-5/6;
b=61/12;
c=-13;
d=107/6;
e=-77/6;
h=15/4;
f_i= e_(x);
f_i_minus_1=e_(x-dx);
f_i_minus_2=e_(x-2*dx);
f_i_minus_3=e_(x-3*dx);
f_i_minus_4=e_(x-4*dx);
f_i_minus_5=e_(x-5*dx);
BDS_fourth_order=((a*f_i_minus_5)+(b*f_i_minus_4)+(c*f_i_minus_3)+(d*f_i_minus_2)+(e*f_i_minus_1)+(h*f_i))/(dx^2);
end
function output = e_(x)
output =exp(x)*cos(x);
end
% slope
function slope= gradient(x1,y1,x2,y2)
slope=(y2-y1)/(x2-x1);
end