[t2, v]=ode45(@ODE4590kg,[0:0.1:300],[0,0]);
xx2_prime = v(:,1);
xx2 = v(:,2);
ma = 90;
gr = 9.81;
we= ma*gr;
co= we/55^2;
dr= co.*(xx2_prime).^2;
bu= 100.*(xx2-8);
re= 3.*(xx2_prime);
if xx2 >= 0 & xx2 <= 8
xx2_double_prime = (we-dr)./ma;
else
xx2_double_prime= (we-dr-bu-re)./ma;
end
figure(2)
subplot(3,2,2)
plot(t2,xx2),title('90 kg ODE 45'),xlabel('Time'),ylabel('Position'),axis([0,300,0,30])
subplot(3,2,4)
plot(t2,xx2_prime),title('90 kg ODE45'),xlabel('Time'),ylabel('Velocity'),axis([0,300,-20,20])
subplot(3,2,6)
plot(t2,xx2_double_prime),title('90 kg ODE 45'),xlabel('Time'),ylabel('Acceleration'),axis([0,300,-20,10])
function [vprime] = ODE4590kg(t2,v )
vprime =[0;0];
m = 90;
g = 9.81;
W= m*g;
c= W/55^2;
D= c.*(v(1)).^2;
B= 100*(v(2)-8);
R= 3.*(v(1));
if v(2) >= 0 && v(2) <= 8
vprime(1) = (W-D)/m;
vprime(2) = v(1);
else
vprime(1) = (W-D-B-R)/m;
vprime(2) = v(1);
end
end
Best Answer