I belive the problem is asking to solve a system of 3 differential equations usi

Question

I belive the problem is asking to solve a system of 3 differential equations using the runge kutta method and simulink method. From what i see the 3 equations are: dy1/dt = y2 * y3 * t ; dy2/dt = -y1 * y3 ; dy3/dt = -0.51 * y1 * y2

The initial conditions appear to be: y1(0) = 0 ; y2(0) = 1 ; y3(0) = 1

Even just the runge kutta method would be helpful.. Thanks!

Solve Example 6.2 on page 111 of textbook, using (1) the 4h order Runge Kutta method and (2) Simulink method, respectively. Please submit your m file and slx file only to Canvas 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)so, y2 (0)=1.0, y3 (0)=1.0 clear; clc; initials(0 . 0 1 . 0 1.0] ; tspan-0.0:0.1:10.0 options odeset ('RelTol',1.0e-6, 'AbsTol, [1.0e-6 t, Y] ode45 (@dydt3,tspan,initial,options) 1.0e-6 1.0e-6]) disp (P) y1=P ( : , 2 ) ; y2-P(:,3) y3-P(:,4) fide fopen ( output . txt , , , w, ) ; fprintf (fid,' y (2) y(3) In') i=1:2:101 fprintf (fid,' for 6.2f %10.2f %10.2f #20.2 n'. end fclose (fid)i plot (t1,Y(:,1) , t 1 , Y ( : , 2 ) , ,-.,,t1,y(:,3),'--'), x1abel('t'), ylabel (Y(1),Y (2),Y(3)),title ('Y 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) Yprime (2)Y(1)Y (3) Yprime (3)-0.51*Y (1)Y (2)

Explanation / Answer

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 x %.3f z %.3f ',t,X,Z);
fprintf('W %.3f ',W);
fprintf('Phi %.3f ',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(' ');

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(' ');

end
end

Related Questions

Navigate

• Browse (All)
• Browse I
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.