Failure (adaptive step): maximum number of trials=4 at step=5 dt=-57.84497129232568

[andreasvarga]

Question❓

I wonder if in the above error message a negative value for dt is admissible.

Note: I am using this software to evaluate the value of a function at a certain time value (via integration of an ODE), which in turn is used to integrate another ODE. The computations are fine with other solvers (e.g., Tsit5()), but fail with IRKGaussLegendre()), The above error message resulted in such a situation.

[mikelehu]

Could you send me the code of the ode-problem and all arguments you use in IRKGL16's call?

[andreasvarga]

Unfortunately, it is a very complicated computation, and I am afraid, it would be difficult to provide an easy to understand example. The computation involves the solution of a periodic differential matrix Riccati equation (PDMRE), where the ODE solver is used first to discretize the differential Riccati equation using a Hamiltonian ODE reformulation for which the symplectic solvers are well suited. A multipoint solution is determined, which serves to compute the solution at each moment, starting from the nearest point as initial condition. Thus the evaluation of a time value of the solution involves the integration of the PDMRE, using a suitable solver (e.g., IRKGaussLegendre()).

The computed solution of the PDMRE is used to determine a stabilizing periodic feedback matrix, which is used to compute the eigenvalues of a suitably constructed monodromy matrix. This computation involves a second integration of a periodic ODE, where each time evaluation of the function involves an integration of an ODE. Just during explaining this process, I realize that such an approach is too complicated to be effective and better I have to look for an alternative approach (e.g., using interpolation techniques).
Anyway, even with an alternative ODE solver, the integration fails when calls to IRKGaussLegendre() are performed for function evaluations. With other solver, I encountered no such failure. The example I have is a linearized satellite model, which has to be stabilized on its orbit.

In the meatime, I tested the interpolation based solution and it is working well, so I will probably switch completely to this solution and avoid two level ODE integration. My question was more oriented to find out if such (crazy) computations could lead to some interference between ODE solvers. Thanks for your time reading this message.

[ChrisRackauckas]

Are you solving the ODE in reverse time?

[andreasvarga]

Yes, the ODE to evaluate the function is solved in reverse time. The ODE at the next level is solved in direct time. Could be this an issue?

[ChrisRackauckas]

For this specific algorithm, maybe. We should try and build an MWE with just those features.

[andreasvarga]

The following code could be such an example:

using OrdinaryDiffEq
using IRKGaussLegendre

function x(t)
    u0 = [1.]
    tspan = (1.,t)
    prob = ODEProblem(Ric1ODE!, u0, tspan)
    #sol = solve(prob, IRKGaussLegendre.IRKGL16(maxtrials=4); adaptive = true, reltol = 1.e-10, abstol = 1.e-10, save_everystep = false)
    sol = solve(prob, Tsit5(); reltol = 1.e-10, abstol = 1.e-10, save_everystep = false)
    return sol(t)
end

function Ric1ODE!(du, u, pars, t) 
    du[1] = 2u[1]-1+u[1]^2
end

f(u, p, t) = x(t)[1]^2
u0 = 0.
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8); 
sol(1)

The result of this computation is 0.37231515461728715 .

However, if you comment out the call to IRKGaussLegendre.IRKGL16 and comment the call to Tsit5, the computation runs forever (at least on my Windows machine).

This is not exacly the same error as in the title, but possibly is related to it via some interation between the two solvers.

[andreasvarga]

After a long time (more than 2 minutes) I got the following result:

0.39093356784066763

which is not correct. With AutoVern9(Rodas5(),nonstifftol = 11/10) I got the result 0.3723151546145423.

[andreasvarga]

With the function

function Ric1ODE!(du, u, pars, t) 
    du[1] = -2u[1]-1+u[1]^2
end

I got the value 3.300615689305358 with Tsit5 but with IRKGaussLegendre.IRKGL16 (and a long waiting time) I got the following output:

┌ ``` Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems). └ @ SciMLBase C:\Users\Andreas.julia\packages\SciMLBase\ynHlA\src\integrator_interface.jl:580 1.0833326941500691e24

[mikelehu]

Sorry for the delay. I think I have already detected what error is in the IRKGL16 integrator. I'm going to make some changes to the code and I'll let you know when IRKGL16 works correctly for the example you sent us.

Thank you.

[ChrisRackauckas]

Reverse time integration was addressed.