Problem with ODE Events Options

11 次查看(过去 30 天)
David
David 2019-9-27
I'm trying to simulate a person making a bungee-jump. To do this I am using the ode45 command.
My problem:
So, the person jumps from an unspecified height, and the rope, L=50m, starts to unfold at that moment. However it does not generate a pulling force (S) until the person actually reaches 50m. Then it decelerates the person until the turning point, where it's starting to accelerate the jumper. When he reaches the 50m mark again, it's the same as before, the force in the rope should be zero. So, I need the ode to make the force in the rope == 0 when it's not fully streched, i.e u(1) > 0, since I've defined the origin where the rope is at L. (u(1) is the y-location in the xy-plane).
I was initially trying to make this happen using if-statements (such as: If u(1)>0) for the input-function, but MatLab was not too happy with me using u(1) as an if-condition. So maybe ODE Event Location is the way to go. But I haven't used it before, and feel a bit confused about the arguments and whatnot.
I would try to create the event function something like this:
function [value,isterminal,direction] = some_handle
value=u(1); %I want this to be equal to 0
isterminal=1; % I'm guessing I would need to stop the integration to make the switch S=0
direction=0; % Any direction
end
But I'm unsure how to write the syntax to connect that function to the values of 'u', and also how I would solve the constant "on/off" switches of the rope-force S.
Here is my code currently, where I've removed my failing attempts at making the ode event location work.
Hopefully this was not too much, and I would really appritiate som help.
(There is also air-resistance --> c*v^2)
g=9.81; mass=75; c=0.1; L=50; k=130;
T=50; %Time
Nstep=500;
tspan=linspace(0,T,Nstep);
u0=[L 0]; %B.conditions
f=@(t,u) uprim(u,t,g,mass,k,c,L);
%options=odeset('Event,?)
[t,u]=ode45(f,tspan,u0); % Here I would need additional arg.-options and output.
function f=uprim(u,t,g,mass,k,c,L)
f=[u(2); g - sign(u(2))*c*u(2)^2/mass - k*u(1)/mass]; % gravity - air-resistance - S
%f2=[u(2); g - sign(u(2))*c*u(2)^2/mass]; For when rope-force S == 0
end
plot(t,u(:,1)) % To see how the location varies with time.

回答(1 个)

Raunak Gupta
Raunak Gupta 2019-10-1
Hi,
I understand you are trying to solve the particular ODE problem using ode45 using ODE Event Location. Also from the dummy event function I can understand that you want to detect the time when ode45 stops in term of u(1) (I am assuming that is displacement). So, the rope force S that is mentioned is dependent on u(1) as it’s a string so I suggest putting value parameter as u(1) only.
Here isterminal condition is critical as if you mention it as 1, the first-time value parameter reaches 0, ode45 will stop at that time moment. direction may be mentioned +1 if the event which signifies the value parameter increasing while going to zero. Same applies for direction value –1. So, Event Location serves as a stopping criterion for ode45 function.
Following Code may help:
g=9.81; mass=75; c=0.1; L=50; k=130;
T=50; %Time
Nstep=5000;
tspan=linspace(0,T,Nstep);
u0=[L 0]; % B.conditions
f=@(t,u) uprim(u,t,g,mass,k,c,L);
options = odeset('Events',@myEventsFcn);
[t,u]=ode45(f,tspan,u0,options);
plot(t,u(:,1)) % To see how the location varies with time.
function [value,isterminal,direction] = myEventsFcn(t,u)
value=u(1); %displacement at a given event is zero
isterminal=1; % Stop the ode45 at this moment
direction=1; % Direction which shows position increasing so that string goes into the relax mode
end
function f=uprim(u,t,g,mass,k,c,L)
f=[u(2); g - sign(u(2))*c*u(2)^2/mass - k*u(1)/mass]; % gravity - air-resistance - S
%f2=[u(2); g - sign(u(2))*c*u(2)^2/mass]; For when rope-force S == 0
end
For more Examples you may refer the following:

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by