How can I solve an ODE that involves a condition on the solution?

8 views (last 30 days)
Nifemi Oluwadairo
Nifemi Oluwadairo on 27 Jul 2020
Commented: Walter Roberson on 4 Sep 2020
I am trying to solve the set of equations below. However, I have the addtional condition that if y(t) < u(t), then y(t) = u(t). How can I account for this condition while solving the equations?
u0=0; udot0 = 0; y0=0; ydot0 = 0; % initial conditions
% u_amp, b, k, F and static_load are all constants
function [zdot] = fun(t,z)
u=z(1); udot=z(2); y=z(3); ydot = z(4);
zdot(1)=udot;
zdot(2)=-u_amp*b^2*cos(c - b*t);
zdot(3)=ydot;
if y > u
zdot(4)= (1/m)*((k*(u - y)) - F*sign(ydot - udot) - static_load);
else
zdot(4)= (1/m)*(-F*sign(ydot - udot) - static_load);
end
zdot=zdot';
end
% this is the condition I am trying to add to my function before solving:
%
% for i = 1:length(y)
% if y(i) < u(i)
% y(i) = u(i);
% end
% end

Answers (2)

Ayush Gupta
Ayush Gupta on 4 Sep 2020
The additional condition for y(t) < u(t), then y(t) = u(t), refer to the following workflow:
function [zdot] = fun(t,z)
u=z(1); udot=z(2); y=z(3); ydot = z(4);
zdot(1)=udot;
zdot(2)=-u_amp*b^2*cos(c - b*t);
zdot(3)=ydot;
if(y(t) <= u(t)) % y(t) is the values of y at time t and similarly u(t) is at time t
% Do whatever you want to do when y(t) <= u(t)
zdot(4)= (1/m)*(-F*sign(ydot - udot) - static_load);
else
% Do whatever you want to do when y(t) > u(t)
zdot(4)= (1/m)*((k*(u - y)) - F*sign(ydot - udot) - static_load);
end
zdot=zdot';
end
  1 Comment
Walter Roberson
Walter Roberson on 4 Sep 2020
u and t are scalars. You cannot index them at t, which is not even an integer, so you cannot use if(y(t) <= u(t))
But anyhow, the approach won't work properly: discontinuity in the derivatives will give invalid answers when used with ode45() or similar.

Sign in to comment.


Walter Roberson
Walter Roberson on 4 Sep 2020
However, I have the addtional condition that if y(t) < u(t), then y(t) = u(t)
You will need to redesign your code.
When you call the MATLAB supplied ode*() routines, the mathematics involved in the Runge-Kutta algorithms requires that your equations be continuously differentiable twice within the boundaries of any one ode*() call.
Any test that is the equivalent of min() or max() such as you have with your y(t) = u(t) when y(t) < u(t), is a place where the function has a discontinuity in differentiation. The mathematics used on the ode*() cannot handle that.
If you are fortunate then the ode*() functions will tell you that you have encountered a singularity. If you are not lucky then the ode*() routines will silently produce the wrong answers.
Unless you can find a way to express your boundaries purely in time, you will need to use event functions to detect when you need to apply the condition, and you must signal to terminate the ode*() call at that point. Then you call ode*() again picking up from the point you left off, except with appropriate adjustments to the boundary conditions.
I recommend reading the ballode example to see how to use event functions.

Community Treasure Hunt

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

Start Hunting!