Question

Solve Example 6.2 on page 111 of textbook, using (1) the 4th order Runge Kutta method and (2) Simulink method, respectively. Please submit your MATLAB code m file and slx file You can use the results from the code attached below (from the textbook by ODE45 solver) as a reference Example 6.2 ode45 ex.m % This program solves a system of 3 differential equations % by using ode45 function % y1 (0)=0, y2 (0)=1.0, y3 (0)=1.0 clear; clc initial [0.0 1.0 1.0] tspan 0.0:0.1:10.0; options-odeset ('RelTol',1.0e-6, 'AbsTol', [1.0e-6 1.0e-6 1.0e-6]) [t, Y] ode45 (@dydt3,tspan, initial,options) disp (P) y1 P(,2) y2-P(:,3) y3 P:,4) fid-fopen (output.txt,'w) fprintf (fid,' y (2) y (3) n') for i=1 :2 : 101 fprintf (fid,' %6.2f %10.2f %10.2f 10.2 \ n' end fclose (fid)i plot (tl, Y(),t1,Y(, 2),'-.',ti,Y:,3),), xlabelt, ylabel (Y(1),Y (2), Y (3)),titleY vs. t grid, text (6.o, -1.2, 'y (1), text (7.7, -0.25, 'y (2), text (4.2,0.85, 'y (3) % dydt3.m % functions for example problem function Yprime-dydt3 (t, Y) Yprime-zeros (3,1) Yprime (1)-Y (2) Y (3) t; yprime (2)--Y(1)*Y (3); Yprime (3)-0.51*Y (1) Y (2)

Answer 1

t=(ti:h:tf);
n=length(t);
if t(n)<tf
t(n+1)=tf;
n=n+1;
end
else t=tspan;
end

while(1)
tend=t(np+1);
hh=t(np+1)-t(np);
while(1)
if hh>h,hh=h;
end
while(1)
if tt+hh>tend,hh=tend-tt;
end
k1=dydt(tt,y(i,:),varargin{:})';
ymid=y(1,:)+k1.*hh./2;
k2=dydt(tt+hh/2,ymid,varargin{:})';
ymid=y(i,:)+k2*hh/2;
k3=dydt(tt+hh/2,ymid,varargin{:});
yend=y(i,:);
k4=dydt(tt+hh,yend,varargin{:})';
phi=(k1+2*(k2+k3)+k4)/6;
y(i+1,:)=y(i,:)+phi*hh;
tt=tt+hh;
i=i+1;
if tt>=tend,break,end
end
np=np+1;
tp(np)=tt;
yp(np,:)=y(i,:);
if tt>=tf,break,end
end
plot(tp,yp)
end

clc;
clear all;

h=1.5;   
x = 0:h:3;   
y = zeros(1,length(x));
y(1) = 5;   
F_xy = @(t,r) 3.*exp(-t)-0.4*r;

for i=1:(length(x)-1)
k_1 = F_xy(x(i),y(i));
k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = F_xy((x(i)+h),(y(i)+k_3*h));

y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  
end

function [t,y] = System(f,n,h)
M = [0,1;-3,-4];
y(1, :) = [2, -4];
f =@(t,y)(M*y')';
t(1) = 0;

for j = 1 : n
k1 = h * f(t(j), y(j,:) );
k2 = h * f( t(j) + h/2, y(j,:) + k1/2 );
k3 = h * f( t(j) + h/2, y(j,:) + k2/2 );
k4 = h * f( t(j) + h, y(j,:) + k3 );
y(j+1,:) = y(j,:) + (k1 + 2*k2 + 2*k3 + k4) / 6;
end

clc; % Clears the screen
clear all;

thete=30;
g= 9.81; %Constant

h=0.01; % step size
tfinal = 3;
N=ceil(tfinal/h); %ceil is round up   

t(1) = 0;% initial condition
Vx(1)=0;%initial accleration
X(1)=0;
Vz(1)=0;
Z(1)=20;%initial velocity
fn = 0.866;
Bz = 0.35;
Phi(1)= 30;
W(1)= 0;

%Bxx = ((Z-Z(1))+5.8*cos(thete)+Bz*(sin(Phi)-sin(thete)));
Bxx = 0.5;
%sliding phase
F_X = @(t,X,Vx,Phi) Vx;
F_Z = @(t,Z,Vz,Phi) Vz;
F_Phi= @(t,Phi,W) W;
F_Vx = @(t,X,Vx,Phi)(fn*(cos(Phi)-0.05*(sin(Phi))));
F_Vz = @(t,Z,Vz,Phi)(fn*(sin(Phi)+0.05*(cos(Phi)))-g);
F_W = @(t,Phi,W)((fn*(Bxx+Bz*0.05))/20000);

for i=1:N; % calculation loop

%update time
t(i+1)=t(i)+h;
%update motions main equation
%rotation phase
k_1 = F_X (t(i) ,X(i) ,Vx(i) ,Phi(i));
L_1 = F_Vx(t(i) ,X(i) ,Vx(i) ,Phi(i));
k_2 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1,Phi(i)+0.5*h*L_1);
L_2 = F_Vx(t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1,Phi(i)+0.5*h*L_1);
k_3 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_2,Vx(i)+0.5*h*L_2,Phi(i)+0.5*h*L_2);
L_3 = F_Vx(t(i)+0.5*h,X(i)+0.5*h*k_2,Vx(i)+0.5*h*L_2,Phi(i)+0.5*h*L_2);
k_4 = F_X (t(i)+h ,X(i)+h *k_3,Vx(i)+ h*L_3,Phi(i)+h*L_3);
L_4 = F_Vx(t(i)+h ,X(i)+h *k_3,Vx(i)+ h*L_3,Phi(i)+h*L_3);

k_11 = F_Z (t(i) ,Z(i) ,Vz(i) ,Phi(i));
L_11 = F_Vz(t(i) ,Z(i) ,Vz(i) ,Phi(i));
k_22 = F_Z (t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11,Phi(i)+0.5*h*L_11);
L_22 = F_Vz(t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11,Phi(i)+0.5*h*L_11);
k_33 = F_Z (t(i)+0.5*h,Z(i)+0.5*h*k_22,Vz(i)+0.5*h*L_22,Phi(i)+0.5*h*L_22);
L_33 = F_Vz(t(i)+0.5*h,Z(i)+0.5*h*k_22,Vz(i)+0.5*h*L_22,Phi(i)+0.5*h*L_22);
k_44 = F_Z (t(i)+ h,Z(i)+ h*k_33,Vz(i)+ h*L_33,Phi(i)+h*L_33);
L_44 = F_Vz(t(i)+ h,Z(i)+ h*k_33,Vz(i)+ h*L_33,Phi(i)+h*L_33);

k_5 = F_Phi (t(i) ,Phi(i) ,W(i));
L_5 = F_W (t(i) ,Phi(i) ,W(i));
k_6 = F_Phi (t(i)+0.5*h,Phi(i)+0.5*h*k_5,W(i)+0.5*h*L_5);
L_6 = F_W (t(i)+0.5*h,Phi(i)+0.5*h*k_5,W(i)+0.5*h*L_5);
k_7 = F_Phi (t(i)+0.5*h,Phi(i)+0.5*h*k_6,W(i)+0.5*h*L_6);
L_7 = F_W (t(i)+0.5*h,Phi(i)+0.5*h*k_6,W(i)+0.5*h*L_6);
k_8 = F_Phi (t(i)+ h,Phi(i)+ h*k_7,W(i)+ h*L_7);
L_8 = F_W (t(i)+ h,Phi(i)+ h*k_7,W(i)+ h*L_7);

X(i+1) = X(i) + (h/6)*(k_1+2*k_2+2*k_3+k_4);
Vx(i+1) = Vx(i) + (h/6)*(L_1+2*L_2+2*L_3+L_4);
Z(i+1) = Z(i) + (h/6)*(k_11+2*k_22+2*k_33+k_44);
Vz(i+1) = Vz(i) + (h/6)*(L_11+2*L_22+2*L_33+L_44);
Phi(i+1)= Phi(i) + (h/6)*(k_5+2*k_6+2*k_7+k_8);
W(i+1) = W(i) + (h/6)*(L_5+2*L_6+2*L_7+L_8);

end

figure
plot(X,Z,'--', 'Linewidth', 1.5, 'color', 'red')
xlabel('time')
ylabel('height')
legend('position')
figure
plot(Phi,W,'--', 'Linewidth', 1.5, 'color', 'blue')
xlabel('time')
ylabel('height')
legend('rotation')
figure
plot(t,Vx,'--', 'Linewidth', 1.5, 'color', 'black')
hold on
plot(t,Vz,'--', 'Linewidth', 1.5, 'color', 'yellow')
hold on
plot(t,W,'--','Linewidth',1.5,'color','green')
xlabel('time')
ylabel('height')
legend('Vx','Vz','W')

fprintf('time %.3f\n x %.3f\n z %.3f\n ',t,X,Z);
fprintf('W %.3f\n',W);
fprintf('Phi %.3f\n',Phi);

function RKSystems(a, b, m, N, alpha, f)
h = (b - a)/N;
t = a;
w = zeros(1, m);

for j = 1:m
w(j) = alpha(j);
end

fprintf('t = %.2f;', t);
for i = 1:m
fprintf(' w(%d) = %.10f;', i, w(i));
end
fprintf('\n');

k = zeros(4, m);
for i = 1:N
for j = 1:m
k(1, j) = h*f{j}(t, w);
end

for j = 1:m
k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));
end

for j = 1:m
k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));
end

for j = 1:m
k(4, j) = h*f{j}(t + h, w + k(3));
end

for j = 1:m
w(j) = w(j) + (k(1, j) + 2*k(2, j) + 2*k(3, j) + k(4, j))/6;
end

t = a + i*h;

fprintf('t = %.2f;', t);
for k = 1:m
fprintf(' w(%d) = %.10f;', k, w(k));
end
fprintf('\n');

end
end