31 views (last 30 days)
Show older comments
Blake Roberts on 29 Jun 2021
-
-
Link
Direct link to this question
https://matlabcentral.mathworks.com/matlabcentral/answers/868403-how-do-i-solve-an-ode-that-is-piecewise-defined-and-oscillates-between-the-two-states
Answered: Bisesh on 19 Apr 2024
Accepted Answer: Jan
Open in MATLAB Online
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
% Advance tstart
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
3 Comments Show 1 older commentHide 1 older comment
Show 1 older commentHide 1 older comment
Jan on 29 Jun 2021
Direct link to this comment
https://matlabcentral.mathworks.com/matlabcentral/answers/868403-how-do-i-solve-an-ode-that-is-piecewise-defined-and-oscillates-between-the-two-states#comment_1611973
Do the options contain an event function to stop the integration?
Blake Roberts on 29 Jun 2021
Direct link to this comment
https://matlabcentral.mathworks.com/matlabcentral/answers/868403-how-do-i-solve-an-ode-that-is-piecewise-defined-and-oscillates-between-the-two-states#comment_1612003
Open in MATLAB Online
Currently my event function has the "isterminal" value set to 0 to continue integration, setting it to 1 seems to stop my results from concatenating within the loop.
%% y = 0 event detection
function [position,isterminal,direction] = y0event(t,y)
position = y(1); % The value we want to detect 0 for
isterminal = 0; % Continue integration (1 will halt it)
direction = 0; % 0 can be approached from +/- direction
end
Jan on 30 Jun 2021
Direct link to this comment
https://matlabcentral.mathworks.com/matlabcentral/answers/868403-how-do-i-solve-an-ode-that-is-piecewise-defined-and-oscillates-between-the-two-states#comment_1612608
See my answer as reaktion to this comment.
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Jan on 30 Jun 2021
Open in MATLAB Online
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
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
More Answers (2)
Bjorn Gustavsson on 30 Jun 2021
Edited: Bjorn Gustavsson on 30 Jun 2021
Open in MATLAB Online
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 Comments Show NoneHide None
Show NoneHide None
Jan on 30 Jun 2021
Direct link to this comment
https://matlabcentral.mathworks.com/matlabcentral/answers/868403-how-do-i-solve-an-ode-that-is-piecewise-defined-and-oscillates-between-the-two-states#comment_1612598
Matlab's ODE integrators are valid, if the function to be integrated is smooth (differentiable). Otherwise the stepsize controller has an undefined behavior: Either it stop the integration with an error or the stepsize is set to such a tiny width, that the rounding errors shadow the discontinuity.
In this case, spring_ext_comp can be evaluated multiple times for one step, while some evaluations happen for y(1)<0 and some for y(1) >= 0. At least in this step the calcluated derivative is random. The effect might look small, but the scientifically correct usage would be to use an event function to stop and restart the integration at y(1)==0. Then using your single spring_ext_comp() would be fine, because the two different models are not mixed.
I've seen too many integrations in scientific publications written without asking an experienced mathematician for numerics. Integrating over discontinuities is a classical fault, which would get obvious, if the sensitivity of the result is examined: Small variations of the inputs or parameters must cause only small deviations of the final values. Otherwise the integration is instable and in consequence not a valid simulation.
Bjorn Gustavsson on 30 Jun 2021
Direct link to this comment
https://matlabcentral.mathworks.com/matlabcentral/answers/868403-how-do-i-solve-an-ode-that-is-piecewise-defined-and-oscillates-between-the-two-states#comment_1612718
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.
Sign in to comment.
Bisesh on 19 Apr 2024
%% 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
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
See Also
Categories
MATLABMathematicsNumerical Integration and Differential EquationsOrdinary Differential Equations
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
