# How do I solve an ODE that is piecewise defined, and oscillates between the two states?

1 view (last 30 days)
Blake Roberts on 29 Jun 2021
Commented: Bjorn Gustavsson on 30 Jun 2021
I have a mass-spring-damper model based on an example in a textbook.
It has sinusoidal input/forcing and a stiffness reduction (alpha) when the displacement is greater than 0.
My current approach
• Execute the ODE solver inside a while loop until the predetermined simulation time ( t_f ) is reached
• Use ODE event detection to trigger an event when the position results cross 0 from either side
• The "if" conditional is supposed to detect which side the spring was compressed or extended prior to reaching 0 and then switch to the other state and continue executing the loop with new initial conditions
• "odecheck" is supposed to record a 1 or a 2 depending on which function of first-order ODE's I call
Issues
• My resulting phase plot does not match up with the results from my text and I think the issue is that my code does not actually switch between the extend and compress functions
• odecheck is a string of zeros except for the very last value
Below is an image of the MSD system as well as the code i am using to accumulate the ODE solver outputs and switch between ODE functions to solve.
I've been working on this for quite a while now and I'm stumped. Any help would be greatly appreciated!
while tout < t_f
[T,sol,te,ye,ie] = ode15s(ode_function,[t_0 t_f],y_0,options);
% SOURCE (L46-52): ballode
% Accumulate output. This could be passed as output arguments.
nt = length(T);
tout = [tout; T(2:nt)];
yout = [yout; sol(2:nt,:)];
teout = [teout; te];
yeout = [yeout; ye];
ieout = [ieout; ie];
% Set new initial conditions
y_0(1) = 0;
y_0(2) = sol(nt,2);
% Change ODE based on which side of y=0 its on
if sol(end,1) > 0
ode_function = @extend;
odecheck(nt) = 1;
elseif sol(end,1) <= 0
ode_function = @compress;
odecheck(nt) = 2;
end
t_0 = T(nt);
end
For reference: I don't think there is an issue with my ODE Functions but just in case here are the two ODEs broken down into a set of first order ODEs for the solver to interperet.
%% ODE Functions
function y_dot = extend(t,y)
global m c k a w
u = 10*sin(w*t);
y_dot(1,1) = y(2);
y_dot(2,1) = (u -c*y(2) - k*a*y(1))/m;
end
function y_dot = compress(t,y)
global m c k w
u = 10*sin(w*t);
y_dot(1,1) = y(2);
y_dot(2,1) = (u -c*y(2) - k*y(1))/m;
end
Jan on 30 Jun 2021
See my answer as reaktion to this comment.

Jan on 30 Jun 2021
The purpose of the event function is to stop the integration, such that the outer loop can restart it with the modified model. If isterminal is 0 in all cases, the event function is sleeping. Should it triggerat y(1)==0? Then:
function [position,isterminal,direction] = y0event(t,y)
position = y(1);
isterminal = 1; % !!!
direction = 0;
end

Bjorn Gustavsson on 30 Jun 2021
Edited: Bjorn Gustavsson on 30 Jun 2021
This modification seems to work "sensibly" OK:
Modify the ODEs into one:
function y_dot = spring_ext_comp(t,y,m,c,k,a,w,Amp)
% SPRING_EXT_COMP -
% Scrapped the globals - a man should have some principles in life - this
% is my...
u = Amp*sin(w*t);
y_dot(1,1) = y(2);
if y(1) > 0
y_dot(2,1) = (u -c*y(2) - k*a*y(1))/m;
else
y_dot(2,1) = (u -c*y(2) - k*y(1))/m;
end
Then I simply run for the full time-of-interest:
T = linspace(0,30,3001);
M = 1;
c = 0.1;
K = 10;
A = 1/8;
W = 10;
Amp = 0; % Just to clearly illustrate that the event-capturing properly modifies the spring-force
[T,Y,te,ye,ie] = ode15s(@(t,y) spring_ext_comp(t,y,M,c,K,A,W,Amp),T,[0.1 12],options);
sol = ode15s(@(t,y) spring_ext_comp(t,y,M,c,K,A,W,Amp),T([1 end]),[0.1 12],options);
clf
plot(T,Y(:,2))
hold on
plot(te,ye,'r*')
plot(sol.xe,sol.ye,'c.','markersize',12)
plot(sol.x,sol.y(2,:),'g')
Which looks OK. Then you can test and check for your parameters - it also looked OK for the positions after the initial decaying natural oscillation have damped out with larger amplitudes on one side when I turned up the amplitude of the drivnig force.
HTH
##### 2 CommentsShowHide 1 older comment
Bjorn Gustavsson on 30 Jun 2021
Well, I learnt something today, something dissapointing, but still. I was assuming that if one made the effort of calculating the exact point of an event then that point would be automatically added to the solution, which would've been sufficient for this case. That was a faulty assumption.

R2020a

### Community Treasure Hunt

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

Start Hunting!

Translated by