ode45 integration issues for some choices of physical parameters,

7 次查看（过去 30 天）
Noob
Noob 2024-10-21
评论： Star Strider 2024-10-22
Hi there!
I have a simple, two-dimensional, ode model that solves for the trajectory of a flat plate falling through a fluid flow. Within my ode model, I have also modeled the aerodynamic lift and drag forces on the plate ('wing'). When I go to solve this ode model using ode45, I notice that for certain choices of density, say, that of paper for the wing, and that of still air for the fluid, ode45 seems to happy to solve my equations, for as long of a time interval as I want to specify. However, when I change the densities to be that of aluminum for the wing, and that of still water for the fluid, ode45 is not happy, and gives me the integration tolerance error message:
Warning: Failure at t=2.240647e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(7.105427e-15) at time t.
In trying to troubleshoot this error, I commented out my new lift and drag force models and instead used very simple lift and drag expressions. However, I still get the above integration tolerance error, which suggests that the issues are not with my fluid force model, and that the issue is somewhere else.
How could I best troubleshoot this issue?
Here are some thoughts of mine:
  1. I'm not aware of any division by zero being performed;
  2. the simulations for paper density and still air density look pretty sensible and pass a number of basic sanity checks;
  3. the simulations with a background moving flow field also look pretty sensible and pass a number of basic sanity checks;
  4. the fluid force model is possibly, perhaps likely, non-differentiable at two corner points; and
  5. when I significantly reduce the fluid density, and keep the wing density to be that of aluminum, then the integration seems fine.
Why would changing densities make my ode model hard to integrate by ode45?
Does fluid density, say, that of still water, cause issues for integration? I would have to guess that this is not the real issue.
Update: I fixed the above issue, but still get integration tolerance issues, when I go to vary other parameters.
So, the above question is truly about varying parameters in general (I think), and less specifically about fluid density or wing density.
Thanks in advance,
  2 个评论
Torsten
Torsten 2024-10-21
编辑：Torsten 2024-10-21
How could we tell what the problem might be without trying your code fed with the successful and unsuccessful material properties ?
Noob
Noob 2024-10-22
编辑：Noob 2024-10-22
Hi Torsten!
I think the issue is the non-smooth behavior of my model, at certain points.
I'm going to try to glue functions together, and approximate my model with a simple function in an epsilon-neighborhood of the singularities. Then, take epsilon as small as I can.
What do you think?
Thanks!

请先登录，再进行评论。

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

采纳的回答

Star Strider
Star Strider 2024-10-21
‘I'm not aware of any division by zero being performed;’
The denominator does not have to be equal to zero, it just has to be sufficiently small to reesult in a singularity. That is simply a property of the representation of floating-point numbers. (I am not certain that your system could be ‘stiff’ with a several orders-of-magnitude difference in the system parameters, however if that is the situaiion, usiing an appropriate stiff solver such as ode15s or ode23s could be appropriate, although it may not prevent the singularity.)
  6 个评论
显示 4更早的评论
Noob
Noob 2024-10-22
Thanks Star Strider!
Will do; at the moment, it seems trickier than I had initially thought:
The singularity takes time to develop.
I know where the singularity is, but I don't know when precisely it occurs.
So my idea of approximating my model with a simple function in an epsilon-neighborhood of the singularity will only work if I can pin down the timing of the singularity.

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2024-10-21
编辑：John D'Errico 2024-10-21
This is an absolutely classic case for a stiff system, where for some sets of system parmeters, you have two or more things happening with significantly different time scales. 1000:1 is probably adequate. And that may not even need to be always the case, but at SOME times, the equations may enter a regime where they become stiff. And so you can see the integration perform adequately, but then suddenly, things go to hell.
A classic example might be a model of a system at night, and then the sun rises, introducing rapid temperature transients. And your problem, of a flat plate falling through a fluid medium is surely one I would fear to be stiff. That is, change the viscosity of the medium, or change the attitude of the plate to where it is no longer acting as a parachute. Turbulence can be nasty of course, so under some conditions, you have vortices form around the edges, etc. I have no idea how complex is your model, but I'd expect it to be stiff.
ODE45, being an adaptive routine, will try to reduce the time step to properly handle the rapid transients your model may experience. However, the common symptom is ODE45 fails, because it can't reduce the time step any further. What did it say?
Warning: Failure at t=2.240647e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(7.105427e-15) at time t.
Absolutely classic for a stiff system. No problems for a while, but then ... it gives up.
https://en.wikipedia.org/wiki/Stiff_equation
(Sorry, but how can you possibly be trying to solve this class of problem, and not even consider if your problem is a stiff one?)
The solution in MATLAB is usually to switch to a solver designed to solve stiff problems. This is also the best way to test IF stiffness is the issue. If a stiff solver fixes it, then it was almost surely stiff. That means you want to use any ode solver with a name that ends in s, so ode15s, ode23s are the ones you would look to immediately. (ODE23T is another, one that can handle moderately stiff problems.)
The very nice thing is, USUALLY this resolves the problem, though sometimes even those solvers are inadequate. And the very nice thing is, for the most part, swapping solvers requires you to do little more than change your call to ODE23S or ODE15S from ODE45. Just change the function name.
  1 个评论
Noob
Noob 2024-10-21
Hi John!
I tried using both stiff solvers, and still get the integration error message, and additionally, this message:
Error using deval (line 139)
Attempting to evaluate the solution outside the interval [0.000000e+00, 1.919039e+00] where it is defined.
I think I should probably work on the model term that I feel is non-smooth at a few points. That's the likely cause of issues; that's my suspicion, I think.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by