MATH4446, Project on Numerical Methods for ODEs This project is on numerical methods for solving initial value problems of the form y′ = f(t, y), a ≤ t ≤ b, y(a) = y0 (1) Note: In order to receive...

1 answer below »
The question is aboutRunge-Kutta method RKF45


MATH4446, Project on Numerical Methods for ODEs This project is on numerical methods for solving initial value problems of the form y′ = f(t, y), a ≤ t ≤ b, y(a) = y0 (1) Note: In order to receive full credit print all your results with double precision >format long and label all your figures and explain your plots and results. 1. In this problem you need to write three MATLAB files, fty.m for f(t,y), rk4.m for the classical one-step fourth-order Runge-Kutta method and rk4main.m which sets the interval [a,b], step-size h, initial conditions y0 and calls rk4 to solve the problem (1). (a) Write an m file rk4.m for the classical fourth-order method. The script rk4.m should start with function y = rk4(t0,y0,h,fty) k1 = feval(fty,t0,y0) k2 = ... k3 = ... k4 = ... y = y0 + .... (b) The script file rk4main.m should define a, b, h, y(1) = y0, and t(i), i = 0, 1, ... such that t(1) = a, T (N) = b. Then, it computes the next iterate calling rk4.m y(i+1) = rk4(t(i),y(i),h,’fty’) plot(t,y) ... other commands to compute errors and print results ... (c) the script fty.m defines f(t, y) and start with function f = fty(t,y) f = expression for f(t,y) (d) Run your program by typing rk4main at the MATLAB command line to solve the problem y′ = 3y + exp(3t), 0 ≤ t ≤ 2, y(0) = 0, (2) 1 with h = 0.4, 0.1, 0.01. Plot the exact solution y(t) = t ∗ exp(3t) and all three numerical solutions on the same plot. Use a different marker for each plot and discuss your results. (e) Compute the final global errors at t = 2 for the previous problem with h = 0.4, 0.1, 0.01 and discuss your results and explain the error behavior. 2. This problem consists of writing a MATLAB program for the adaptive Runge-Kutta method RKF45 given in pages 296-298 of the textbook with automatic step-size selec- tion for solving systems of ordinary differential equations of the form dY dt = F (t, Y ), a ≤ t ≤ b, Y (a) = Y0. (3) (a) Write a MATLAB program rkf45.m for adaptive Runge-Kutta-Fehlberg algorithm described in pages 296-298 of the textbook. The function F (t, Y ) should be defined in the matlab file ftys.m. Initialize hmax = 0.2, hmin = 10−5 and in Step 4 of the algorithm w̃ and w are vectors thus |w̃ − w| should be the infinity norm which can be computed using the MATLAB command ’norm(vector,inf)’. Print the number of successful and rejected steps for all problems. (b) Use your program to solve the second-order problem for the swinging pendulum with the angle u(t) versus time t and gravitational acceleration g = −32.17ft/s2 and length l. Print the number of rejected and successful steps. u′′ − g l sin(u) = 0, 0 ≤ t ≤ 10, u(0) = 1.2, u′(0) = 0. (4) i. Transform the equation of (4) into a system of first order differential equations and define F (t, y) in ftys function f = ftys(t,y) f(1) = ... f(2) = ... Note that y0 is a vector, feval in rkf45.m returns a vector. Thus, for each t(i) you have a vector [y(i, 1), y(i, 2)] approximating [y1(t(i)), y2(ti))]. ii. Solve (4) with l = 0.2 using TOL = 0.1 and TOL = 0.0001. Plot the two numerical solutions of u(t) on the same plot. Create another figure and plot the two approximations of u′(t). iii. Discuss your numerical results. (c) Solve problem 3.c, page 338, with TOL = 0.0001 and plot the true solution y(t) and its numerical approximation. Print the number of rejected steps and the number of successful steps. 2
Answered Same DayMay 04, 2021

Answer To: MATH4446, Project on Numerical Methods for ODEs This project is on numerical methods for solving...

Rahul answered on May 08 2021
152 Votes
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 =...
SOLUTION.PDF

Answer To This Question Is Available To Download

Related Questions & Answers

More Questions »

Submit New Assignment

Copy and Paste Your Assignment Here