Ode45, evolution giving NaNs after some time, with no evidence of singularities in the equations

Hello,

Attached is a simple integrator for four coupled first order ode s.

after a certain time, the 4-vector solution turns into NaNs with no evidence of singularities in the equations.

This happens after a certain range where the solutions (starting as real) are written in complex notation although the imaginary part is zero. (why?)

can anyone please have a look and let me know if there's something im missing?

[darova]

I tried to plot imiginary

plot(t1,imag(x(:,2)))

Result

[Mohammad Farhat]

yes I said the notation at somepoint (if you check the entries of the solutions), for no apparent reason to me at least, becomes a complex notation, although as you've seen the solution is still purely real.

Answer 1

[Walter Roberson]

Your first order, e1, gets to be more than 1.0, and then the sqrt(1-x^2) gets to be complex. Problems accumulate from there.

What are the original equations that you are trying to deal with? Do you have the Symbolic Toolbox ?

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/274160-when-use-ode45-there-nan-output-why-and-how-to-solve-it#comment_823368

[Mohammad Farhat]

I see. When I looked at the solution, the last values before the nans were were pure real (as darova plotted above, although written in complex notation but with 0i) so I assumed e is not going beyond 1.0 ,

The original equations I'm trying to solve are the four Eqs. 24-27 giving the evolution figure next to them in figure 1 in the following paper:

https://arxiv.org/pdf/1409.6717.pdf

Do you have any advice on how to fix this?

[Walter Roberson]

I did a bunch of tracing, and figured out what is going on.

ode45 is an adaptive-step system. It moves forward and tends to adjust step size upwards until it encounters a location where the error would be more than the configured tolerance, and then it backs up on the step.

Sometimes during this process, it projects an e value larger than 1. That e value leads to complex-valued results. If this occurs during the middle of the adaptive sequence, the complex results can end up getting fed in, and that leads to even larger complex results. Most of the time, MATLAB would then simply detect that the error was too large and would then back up and try again.

However, especially as O gets more and more negative, if the 4th of 6 function evaluations in an ode45 cycle turns out "a little bit" complex, the trial values for the 6th function evaluation can contain values such that the sin() or cos() explode (because sin and cos of complex numbers can be defined in terms of sinh() and cosh() and those involve exp() of the individual real and imaginary coefficients.) That leads to complex infinities, and those quickly lead to nan.

There is a work-around: at the end of your rhs routine, add

    if any(isnan(dxdt))
        dxdt(:) = 1e10;
    end

This replaces the nan with large values, large enough to trigger the error correction.

It takes a while to get through the entire time span you have specified.

Because of this correction process, even though individual evaluations can end up being complex valued, those all happen to get rejected as generating too abrupt of a change, so you eventually end up with final values that are all real-valued. But you do need to work around the temporary nan values triggered by the step being too large at the time.

[Mohammad Farhat]

Thanks a lot for this explanation.

So this removed the nans (but still I dont have the needed results as in the figure..)

What do you think about using an event funtion for e ( x(1) ) so that it doesnt exceed 1? (this is a physical constraint )

[darova]

Thank you. This is something new

[Walter Roberson]

I do not think an event function is appropriate here. Projection and failed steps is normal for ode45() .

Setting a MaxStep in the options might potentially help.