Answer To: MATH4446, Project on Numerical Methods for ODEs This project is on numerical methods for solving...
Rahul answered on May 08 2021
math_050520/fty.m
% Matlab function to change string to function in this independent variable
% is x
% Input f in string format
% Ouput function_handle
function fty = fty()
%f = '(5*(t.^2) - y)/(exp(t+y))'
f = '3*y + exp(3*t)';
str = strcat('@','(t,y)',f);
fty = str2func(str);
end
math_050520/ftys.m
function [f1,f2,y0] = ftys(t,y)
g = -32.17;
l = 0.2;
% syms u(t)
% eqn = diff(u,2) - (g/l)*sin(u) ==0;
% [eq_f,vars] = reduceDifferentialOrder(eqn,u(t))
% fun_1 = eq_f(1);
% func_2 = eq_f(2);
% a = matlabFunction(fun_1)
% a_s = func2str(a)
% a_s = a_s(1:(length(a_s)-15));
% str = strrep( a_s, 'u(t)', 'y' );
% str = strrep(str,'(t)','(t,y,u1)(-1)*');
% [V,s] = odeToVectorField(eqn);
% M = matlabFunction(V,'vars',{'t','Y'}
f = @(t,u,u1)(165.8)*sin(u);
f2 = f
f1 = @(t,u,u1)u1;
y0 = [1.2;0];
end
math_050520/h_0.01.png
math_050520/h_0.1.png
math_050520/h_0.4.png
math_050520/part_1_solve.m
% Matlab Code to solve the function by Ranga Kutta numerical integration
% method
% Intial define the variable t and y
format long e % as per instruction to get double precison result
% f = '3*y + exp(3*t)'; Given function
fty = fty();% Change the string function to function handle
a = 0;% Lower limit of t
b = 2;% Upper Limit of t
h = [0.4,0.1,0.01];% define h matrix so that easy to find the solution at different h value
h_1 = 0.4;
h_2 = 0.1;
h_3 = 0.01
% h = 0.4
% Analysis the relation between numerical solution and exact solution for h = 0.4
y0 = 0;% Initial y value
% n = 1:length(h) % Using the Loop to find solution at different h value
t_1 = a:h_1:b;% divide the domain into steps
y_1 = rk4main(a,b,h_1,y0,fty);% find the value of y
figure(1);
a_1_1 = zeros(length(t_1),1);
for i = 1:length(t_1);% Exact Solution
a_1_1(i,1) =(t_1(i)*exp(3*t_1(i)));
end
error_1 = a_1_1(end)-y_1(end)
yyaxis left
plot(t_1',y_1,'b--o');
ylabel('y');
yyaxis right
plot(t_1',a_1_1,'r--o');
title('When h Value is 0.4');
ylabel('y');
xlabel('t');
legend('y approx','y exact')
% Analysis the relation between numerical solution and exact solution for h = 0.1
% y0 = 0;% Initial y value
% n = 1:length(h) % Using the Loop to find solution at different h value
t_2 = a:h_2:b;% divide the domain into steps
y_2 = rk4main(a,b,h_2,y0,fty);% find the value of y
figure(2);
a_1_2 = zeros(length(t_2),1);
for i = 1:length(t_2);% Exact Solution
a_1_2(i,1) =(t_2(i)*exp(3*t_2(i)));
end
error_2 = a_1_2(end)-y_2(end)
yyaxis left
plot(t_2',y_2,'ro');
ylabel('y');
yyaxis right
plot(t_2',a_1_2,'k--*');
title('When h Value is 0.1');
ylabel('y');
xlabel('t');
legend('y approx','y exact')
% Analysis the relation between numerical solution and exact solution for h = 0.01
% y0 = 0;% Initial y value
% n = 1:length(h) % Using the Loop to find solution at different h value
t_3 = a:h_3:b;% divide the domain into steps
y_3 = rk4main(a,b,h_3,y0,fty);% find the value of y
figure(3);
a_1_3 = zeros(length(t_3),1);
for i = 1:length(t_3);% Exact Solution
a_1_3(i,1) =(t_3(i)*exp(3*t_3(i)));
end
error_3 = a_1_3(end)-y_3(end)
yyaxis left
plot(t_3',y_3,'r--o');
ylabel('y');
yyaxis right
plot(t_3',a_1_3,'k--*');
title('When h Value is 0.01');
ylabel('y');
xlabel('t');
legend('y approx','y exact')
disp(['The absolute global error when h is 0.4 = ',num2str(error_1)]);
disp(['The absolute global error when h is 0.1 = ',num2str(error_2)]);
disp(['The absolute global error when h is 0.01 = ',num2str(error_3)]);
math_050520/part_2a.m
% Matlab Code to solve the 296-298 example
y = 'y - (t^2) +1'
ftys = @(t,y)y - (t^2) +1
a = 0;
b = 2;
alpha = 0.5;
TOL = 10e-5;
w = RKF45(a,b,alpha,TOL,ftys);
math_050520/part_3c.asv
% Matlab Code to solve Runga Kutta fehleberg
% Matlab Code is written to find the numerical integration by Runga Kutta
% Fehlberg
f1 = @(t,y,y1,y2)y1
f2 = @(t,y,y1,y2)y2
f3 = @(t,y,y1,y2)(-2*y2 + y1 + 2*y + exp(t))
y0_1 = 1;
y0_2 = 2;
y0_3 = 0;
a = 0; % Lower limit of t
b = 3;% Upper Limit of t
alpha1 = y0_1;
alpha2 = y0_2;% Initial y value
alpha3 = y0_3;
TOL = 0.0001; % Tolerance
h = 0.2;
t(1) = a;% staring t value
w1(1) = alpha1;% startng y value
w2(1) = alpha2;
w3(1) = alpha3;
FLAG = 1;
% ftys = @(t,y)3*y + exp(3*t);
fprintf(' t w\n');
ss = 0;
i = 1;
rs = 0;
fprintf('%5.4f %11.8f \n', t, w1);
while FLAG ==1
K1 = h*feval(f1,t(i),w1(i),w2(i),w3(i));
L1 = h*feval(f2,t(i),w1(i),w2(i),w3(i));
M1 = h*feval(f3,t(i),w1(i),w2(i),w3(i));
K2 = h*feval(f1,t(i) + (0.25*h),w1(i) + (0.25*K1),w2(i) + (0.25*L1),w3(i) + (0.25*M1));
L2 = h*feval(f2,t(i) + (0.25*h),w1(i) + (0.25*K1),w2(i) + (0.25*L1),w3(i) + (0.25*M1));
M2 = h*feval(f3,t(i) + (0.25*h),w1(i) + (0.25*K1),w2(i) + (0.25*L1),w3(i) + (0.25*M1));
K3 = h*feval(f1,t(i) + (3*h/8),w1(i) + (3*K1/32) + (9*K2/32),w2(i) + (3*L1/32) + (9*L2/32),w3(i) + (3*M1/32) + (9*M2/32));
L3 = h*feval(f2,t(i) + (3*h/8),w1(i) + (3*K1/32) + (9*K2/32),w2(i) + (3*L1/32) + (9*L2/32),w3(i) + (3*M1/32) + (9*M2/32));
M3 = h*feval(f3,t(i) + (3*h/8),w1(i) + (3*K1/32) + (9*K2/32),w2(i) + (3*L1/32) + (9*L2/32),w3(i) + (3*M1/32) + (9*M2/32));
K4 = h*feval(f1,t(i) + (12*h/13),w1(i) + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197),w2(i) + (1932*L1/2197) - (7200*L2/2197) + (7296*L3/2197),w3(i) + (1932*M1/2197) - (7200*M2/2197) + (7296*M3/2197));
L4 = h*feval(f2,t(i) + (12*h/13),w1(i) + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197),w2(i) + (1932*L1/2197) - (7200*L2/2197) + (7296*L3/2197),w3(i) + (1932*M1/2197) - (7200*M2/2197) + (7296*M3/2197));
M4 = h*feval(f3,t(i) + (12*h/13),w1(i) + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197),w2(i) + (1932*L1/2197) - (7200*L2/2197) + (7296*L3/2197),w3(i) + (1932*M1/2197) - (7200*M2/2197) + (7296*M3/2197));
K5 = h*feval(f1,t(i) + h,w1(i) + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104),w2(i) + (439*L1/216) - (8*L2) + (3680*L3/513) - (845*L4/4104),w3(i) + (439*M1/216) - (8*M2) + (3680*M3/513) - (845*M4/4104));
L5 = h*feval(f2,t(i) + h,w1(i) + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104),w2(i) + (439*L1/216) - (8*L2) + (3680*L3/513) - (845*L4/4104),w3(i) + (439*M1/216) - (8*M2) + (3680*M3/513) - (845*M4/4104));
M5 = h*feval(f3,t(i) + h,w1(i) + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104),w2(i) + (439*L1/216) - (8*L2) + (3680*L3/513) - (845*L4/4104),w3(i) + (439*M1/216) - (8*M2) + (3680*M3/513) - (845*M4/4104));
K6 = h*feval(f1,t(i) + (h/2),w1(i) - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40),w2(i) - (8*L1/27) + (2*L2) - (3544*L3/2565) + (1859*L4/4104) - (11*L5/40),w3(i) - (8*M1/27) + (2*M2) - (3544*M3/2565) + (1859*M4/4104) - (11*M5/40));
L6 = h*feval(f2,t(i) + (h/2),w1(i) - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40),w2(i) - (8*L1/27) + (2*L2) - (3544*L3/2565) + (1859*L4/4104) - (11*L5/40),w3(i) - (8*M1/27) + (2*M2) - (3544*M3/2565) + (1859*M4/4104) - (11*M5/40));
K6 = h*feval(f3,t(i) + (h/2),w1(i) - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40),w2(i) - (8*L1/27) + (2*L2) - (3544*L3/2565) + (1859*L4/4104) - (11*L5/40),w3(i) - (8*M1/27) + (2*M2) - (3544*M3/2565) + (1859*M4/4104) - (11*M5/40));
wc = w1(i) + (16*K1/135) + (6656*K3/12825) + (28561*K4/56430) - (9*K5/50) + (2*K6/55);
we = w1(i) + (25*K1/216) + (1408*K3/2565) + (2197*K4/4104) - (K5/5);
vect = wc - we;
nor = norm(vect,inf);
R = nor/h;
% R = (1/h)*abs((K1/360) - (128*K3/4275) - (2197*K4/75240) + (1*K5/50) + (2*K6/55))
if R <= TOL
t(i+1) = t(i) + h;
w1(i+1) = w1(i) + (25*K1/216) + (1408*K3/2565) + (2197*K4/4104) - (K5/5);
w2(i+1) = w2(i) + (25*L1/216) + (1408*L3/2565) + (2197*L4/4104) - (L5/5);
w3(i+1) = w3(i) + (25*M1/216) + (1408*M3/2565) + (2197*M4/4104) - (M5/5);
ss = ss + 1;
else
t(i+1) = t(i);
w1(i+1) = w1(i);
w2(i+1) = w2(i);
w3(i+1) = w3(i)
end
delta = (0.84)*((TOL/R)^0.25);
if delta <= 0.1
h = 0.1*h;
elseif delta >= 4
h = 4*h;
else
h = delta*h;
end
if t(i)>=b
FLAG = 0;
end
fprintf('%5.4f %11.8f \n', t(i+1), w1(i+1));
% plot(t,w1);
% hold on
i = i + 1;
end
disp(['Number of Rejected Steps = ',num2str(rs)]);
disp(['Number of Successfull Steps = ',num2str(ss)]);
for i_1 = 1: length(t)
a_1(i_1) = ((43/36)*exp(t(i_1))) + (0.25*exp(-t(i_1))) - ((4/9)*exp(-2*t(i_1))) + ((1/6)*(t(i_1)*(exp(t(i_1)))));
end
plot(t,w1,'r',t,a_1,'b');
legend('approximat Solution','Exact Solution');
math_050520/part_3c.m
% Matlab Code to solve Runga Kutta fehleberg
% Matlab Code is written to find the numerical integration by Runga Kutta
% Fehlberg
f1 = @(t,y,y1,y2)y1;
f2 = @(t,y,y1,y2)y2;
f3 = @(t,y,y1,y2)(-2*y2 + y1 + 2*y + exp(t));
y0_1 = 1;
y0_2 = 2;
y0_3 = 0;
a = 0; % Lower limit of t
b = 3;% Upper Limit of t
alpha1 = y0_1;
alpha2 = y0_2;% Initial y value
alpha3 = y0_3;
TOL = 0.0001; % Tolerance
h = 0.2;
t(1) = a;% staring t value
w1(1) = alpha1;% startng y value
w2(1) = alpha2;
w3(1) = alpha3;
FLAG = 1;
% ftys = @(t,y)3*y + exp(3*t);
fprintf(' t w\n');
ss = 0;
i = 1;
rs = 0;
fprintf('%5.4f %11.8f \n', t, w1);
while FLAG ==1
K1 = h*feval(f1,t(i),w1(i),w2(i),w3(i));
L1 = h*feval(f2,t(i),w1(i),w2(i),w3(i));
M1 = h*feval(f3,t(i),w1(i),w2(i),w3(i));
K2 = h*feval(f1,t(i) + (0.25*h),w1(i) + (0.25*K1),w2(i) + (0.25*L1),w3(i) + (0.25*M1));
L2 = h*feval(f2,t(i) + (0.25*h),w1(i) + (0.25*K1),w2(i) + (0.25*L1),w3(i) + (0.25*M1));
M2 = h*feval(f3,t(i) + (0.25*h),w1(i) + (0.25*K1),w2(i) + (0.25*L1),w3(i) + (0.25*M1));
K3 = h*feval(f1,t(i) + (3*h/8),w1(i) + (3*K1/32) + (9*K2/32),w2(i) + (3*L1/32) + (9*L2/32),w3(i) + (3*M1/32) + (9*M2/32));
L3 = h*feval(f2,t(i) + (3*h/8),w1(i) + (3*K1/32) + (9*K2/32),w2(i) + (3*L1/32) + (9*L2/32),w3(i) + (3*M1/32) + (9*M2/32));
M3 = h*feval(f3,t(i) + (3*h/8),w1(i) + (3*K1/32) + (9*K2/32),w2(i) + (3*L1/32) + (9*L2/32),w3(i) + (3*M1/32) + (9*M2/32));
K4 = h*feval(f1,t(i) + (12*h/13),w1(i) + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197),w2(i) + (1932*L1/2197) - (7200*L2/2197) + (7296*L3/2197),w3(i) + (1932*M1/2197) - (7200*M2/2197) + (7296*M3/2197));
L4 = h*feval(f2,t(i) + (12*h/13),w1(i) + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197),w2(i) + (1932*L1/2197) - (7200*L2/2197) + (7296*L3/2197),w3(i) + (1932*M1/2197) - (7200*M2/2197) + (7296*M3/2197));
M4 = h*feval(f3,t(i) + (12*h/13),w1(i) + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197),w2(i) + (1932*L1/2197) - (7200*L2/2197) + (7296*L3/2197),w3(i) + (1932*M1/2197) - (7200*M2/2197) + (7296*M3/2197));
K5 = h*feval(f1,t(i) + h,w1(i) + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104),w2(i) + (439*L1/216) - (8*L2) + (3680*L3/513) - (845*L4/4104),w3(i) + (439*M1/216) - (8*M2) + (3680*M3/513) - (845*M4/4104));
L5 = h*feval(f2,t(i) + h,w1(i) + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104),w2(i) + (439*L1/216) - (8*L2) + (3680*L3/513) - (845*L4/4104),w3(i) + (439*M1/216) - (8*M2) + (3680*M3/513) - (845*M4/4104));
M5 = h*feval(f3,t(i) + h,w1(i) + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104),w2(i) + (439*L1/216) - (8*L2) + (3680*L3/513) - (845*L4/4104),w3(i) + (439*M1/216) - (8*M2) + (3680*M3/513) - (845*M4/4104));
K6 = h*feval(f1,t(i) + (h/2),w1(i) - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40),w2(i) - (8*L1/27) + (2*L2) - (3544*L3/2565) + (1859*L4/4104) - (11*L5/40),w3(i) - (8*M1/27) + (2*M2) - (3544*M3/2565) + (1859*M4/4104) - (11*M5/40));
L6 = h*feval(f2,t(i) + (h/2),w1(i) - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40),w2(i) - (8*L1/27) + (2*L2) - (3544*L3/2565) + (1859*L4/4104) - (11*L5/40),w3(i) - (8*M1/27) + (2*M2) - (3544*M3/2565) + (1859*M4/4104) - (11*M5/40));
K6 = h*feval(f3,t(i) + (h/2),w1(i) - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40),w2(i) - (8*L1/27) + (2*L2) - (3544*L3/2565) + (1859*L4/4104) - (11*L5/40),w3(i) - (8*M1/27) + (2*M2) - (3544*M3/2565) + (1859*M4/4104) - (11*M5/40));
wc = w1(i) + (16*K1/135) + (6656*K3/12825) + (28561*K4/56430) - (9*K5/50) + (2*K6/55);
we = w1(i) + (25*K1/216) + (1408*K3/2565) + (2197*K4/4104) - (K5/5);
vect = wc - we;
nor = norm(vect,inf);
R = nor/h;
% R = (1/h)*abs((K1/360) - (128*K3/4275) - (2197*K4/75240) + (1*K5/50) + (2*K6/55))
% % if R <= TOL
t(i+1) = t(i) + h;
w1(i+1) = w1(i) + (25*K1/216) + (1408*K3/2565) + (2197*K4/4104) - (K5/5);
w2(i+1) = w2(i) + (25*L1/216) + (1408*L3/2565) + (2197*L4/4104) - (L5/5);
w3(i+1) = w3(i) + (25*M1/216) + (1408*M3/2565) + (2197*M4/4104) - (M5/5);
ss = ss + 1;
% else
% t(i+1) = t(i);
% w1(i+1) = w1(i);
% w2(i+1) = w2(i);
% w3(i+1) = w3(i);
% end
% delta = (0.84)*((TOL/R)^0.25);
% if delta <= 0.1
% h = 0.1*h;
% elseif delta >= 4
% h = 4*h;
% else
% h = delta*h;
% end
if t(i)>=b
FLAG = 0;
end
fprintf('%5.4f %11.8f \n', t(i+1), w1(i+1));
% plot(t,w1);
% hold on
i = i + 1;
end
disp(['Number of Rejected Steps = ',num2str(rs)]);
disp(['Number of Successfull Steps = ',num2str(ss)]);
for i_1 = 1: length(t)
a_1(i_1) = ((43/36)*exp(t(i_1))) + (0.25*exp(-t(i_1))) - ((4/9)*exp(-2*t(i_1))) + ((1/6)*(t(i_1)*(exp(t(i_1)))));
end
plot(t,w1,'r',t,a_1,'k--o');
legend('approximat Solution','Exact Solution');
math_050520/part_c.png
math_050520/rk4.m
%Matlab Code to solve the ranga kutta fourth order equation
% Input:
% y0 = Initial y value
% t0 = Initial t value
% h = step value;
% fty = function handle which comes from the fty.m
function y = rk4(t0,y0,h,fty)
%t0 = t(n-1);y0 = y(n-1,1);h = 1.5;
k1 = h*feval(fty,t0,y0);
k2 = h*feval(fty,t0 + (0.5*h),y0 + (0.5*k1));
k3 = h*feval(fty,t0 + (0.5*h),y0 + (0.5*k2));
k4 = h*feval(fty,t0 + (h),y0 + (k3));
y = y0 + ((k1 + 2*k2 + 2*k3 + k4)/6); % Valye of y after one step
end
% Output: y value after one step
math_050520/rk4main.m
% Matlab Code for the Ranga kutta numerical integration by doing lots of steps
% Input:
% a = Lower limit t value
% b = Upper limit t value
% h = step value;
% y0 = initial y valye
function y = rk4main(a,b,h,y0,fty)
t = a:h:b;%make a matrix of t value after dividing into number of steps
n = length(t);% Number of interval
y = zeros(n,1);% Initila y value consideration
y(1,1) = y0;% y value at initial steps
for i = 2:n;% find the value of y at all steps
y(i,1) = rk4(t(i-1),y(i-1,1),h,fty);
end
plot(t,y);% Plot between t and y
end
% Outcome:
% y value at given number of steps
math_050520/RKF45.m
% Matlab Code is written to find the numerical integration by Runga Kutta
% Fehlberg
function w = RKF45(a,b,alpha,TOL,ftys)
syms t y
% Input
% ftys = ftys(t,y);% Input as function
% a = 0; % Lower limit of t
% b = 2;% Upper Limit of t
% alpha = 0.25;% Initial y value
% TOL = 1e-5; % Tolerance
hmax = 0.25; % maiximum step
hmin = 0.01;% Minimum step
t = a;% staring t value
w = alpha;% startng y value
h = hmax;
FLAG = 1;
% ftys = @(t,y)3*y + exp(3*t);
fprintf(' t w\n');
ss = 0;
rs = 0;
fprintf('%5.4f %11.8f \n', t, w);
while FLAG ==1
K1 = h*feval(ftys,t,w);
K2 = h*feval(ftys,t + (0.25*h),w + (0.25*K1));
K3 = h*feval(ftys,t + (3*h/8),w + (3*K1/32) + (9*K2/32));
K4 = h*feval(ftys,t + (12*h/13),w + (1932*K1/2197) - (7200*K2/2197) + (7296*K3/2197));
K5 = h*feval(ftys,t + h,w + (439*K1/216) - (8*K2) + (3680*K3/513) - (845*K4/4104));
K6 = h*feval(ftys,t + (h/2),w - (8*K1/27) + (2*K2) - (3544*K3/2565) + (1859*K4/4104) - (11*K5/40));
wc = w + (16*K1/135) + (6656*K3/12825) + (28561*K4/56430) - (9*K5/50) + (2*K6/55);
we = w + (25*K1/216) + (1408*K3/2565) + (2197*K4/4104) - (K5/5);
vect = wc - we;
nor = norm(vect,inf);
R = nor/h;
% R = (1/h)*abs((K1/360) - (128*K3/4275) - (2197*K4/75240) + (1*K5/50) + (2*K6/55))
if R <= TOL
t = t + h;
w = w + (25*K1/216) + (1408*K3/2565) + (2197*K4/4104) - (K5/5);
ss = ss + 1;
else
rs = rs +1;
end
delta = (0.84)*((TOL/R)^0.25);
if delta <= 0.1
h = 0.1*h;
elseif delta >= 4
h = 4*h;
else
h = delta*h;
end
if (h > hmax)
h = hmax;
end
if t>=b
FLAG = 0;
elseif (t + h) > b
h = b - t;
elseif h < hmin
FLAG = 0;
end
fprintf('%5.4f %11.8f \n', t, w);
end
disp(['Number of Rejected Steps = ',num2str(rs)]);
disp(['Number of Successfull Steps = ',num2str(ss)]);
math_050520/Solution.pdf
Study Of Numerical Integration
1. Matlab File for fty.m
% Matlab function to change string to function in this independent variable
% is x
% Input f in string format
% Ouput function_handle
function fty = fty()
%f = '(5*(t.^2) - y)/(exp(t+y))'
f = '3*y + exp(3*t)';
str = strcat('@','(t,y)',f);
fty =...