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

8 次查看（过去 30 天）
Nifemi Oluwadairo
Nifemi Oluwadairo 2020-7-27
评论： Walter Roberson 2020-9-4
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

请先登录，再回答此问题。

回答（2 个）

Ayush Gupta
Ayush Gupta 2020-9-4
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 个评论
Walter Roberson
Walter Roberson 2020-9-4
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2020-9-4
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.

类别

在 Help Center 和 File 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