The question is aboutRunge-Kutta method RKF45MATH4446, Project on Numerical Methods for ODEs This...

Question

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

project4-0efca0zj.pdf 338-1mlf3r15-fh45fnyd.pdf 296-298-rgzrflq5-vokxw2re.pdf

Rahul · Accepted Answer

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
');
ss = 0;
i = 1;
rs = 0;
fprintf('%5.4f  %11.8f  
', 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,

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...

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

Answer To This Question Is Available To Download

Related Questions & Answers

Submit New Assignment