ode45 integration issues for some choices of physical parameters,

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,

[Torsten]

How could we tell what the problem might be without trying your code fed with the successful and unsuccessful material properties ?

[Noob]

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!

Answer 1

[Star Strider]

‘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.)

[Noob]

Hi there!

I just tried both stiff solvers, and still get the integration issue, along with this new, additional error message:

Error using deval (line 139)

Attempting to evaluate the solution outside the interval [0.000000e+00, 1.919039e+00] where it is defined.

Does that mean my model is non-smooth, or not differentiable, or not twice-differentiable (C^2 function), for at least one point?

In my making a bunch of figures before I wrote the ode model, I am aware that a term within my ode model is non-smooth: its curve has two kinks / corner points on it, which I have been unable to resolve, so far.

So, perhaps my equations aren't stiff, so much as my equations are non-smooth (because a model term is non-smooth)?

I just read about stiff equations -- about stiffness being about computational inefficiency. Stiffness does not appear to be about the model's smoothness.

Thanks!

[Star Strider]

My pleasure!

If it is has a significant discontinuity (is non-differentiable) at any point, the usual approach is to integrate up to the first discontinuity, stop the solver, re-start the solver using the final values for the variables and time as the initial values for the next step, continue to the second discontinuity, and repeat that process with the new final values for the last part of the integration. If you know that there are three separate parts of the integration, you can use a for loop to start and stop the integration at the appropriate times.

Use the tspan vector to define the times delineating the segments if you know what the times are that define them. Otherwise, you may need to use an Events function. See ODE Event Location for those details.

[Noob]

How about I make my model smooth, or as smooth as possible, by gluing functions together, and approximating my model in an epsilon-neighborhood of the singularity with a simple function? That is what I think I'll try next?

[Star Strider]

I know nothing about your model. If you believe that approach to be a solution to the problem you are curreently having, do it!

Please report back here on the result.

[Noob]

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.

[Star Strider]

My pleasure!

That, or could experiment with an event.

Answer 2

[John D'Errico]

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.

[Noob]

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.