1
$\begingroup$

I am trying to numerically integrate a differential equation using scipy.integrate.RK45 and/or scipy.integrate.LSODA.

Now, I am trying to fix the integration step sizes of both solvers. This, however, does not work out. I have tried the following:

1.) Setting the "h_abs" method to my desired step size BEFORE every call of "step()". This however gets ignored by the solver.

2.) I have tried to set my "rtol" and "atol" values very high. This too gets ignored.

3.) Both 1.) and 2.) simultaneously.

I have also noticed a strange behaviour of my solvers. The problem at hand is a PDE to which I apply a finite difference scheme in the spatial domain such that it becomes an ODE in time. For a large courant factor, the solvers accept my fixed step sizes with the above methods. When my courant factor becomes smaller, the solvers start to randomly adjust their step sizes. This is not only weird but also very annoying, since I cannot really work with some solver that does what it wants without control.

Any idea what causes this behaviour and how I can circumvent it?

$\endgroup$
4
  • $\begingroup$ Turning a PDE into an ODE system via a fixed space grid is called "method of lines". $\endgroup$ Nov 26, 2022 at 11:46
  • $\begingroup$ Any rules you apply deal with the inherent stiffness of the PDE, due to the condition number of the matrix of the discretized space derivatives going to infinity with decreasing step size of the space grid. With variable-step methods, the step size controller "feels" the boundaries of the stability region via error estimates, so a sub-critical step size is enforced automatically. So the problem you are trying to solve may not actually exist. $\endgroup$ Nov 26, 2022 at 11:53
  • $\begingroup$ Just to avoid any misconception: Do I understand you right, that by applying the FDA scheme my ODE becomes stiff when decreasing the step size? This would explain a lot; I tried to implement an RK4 solver myself such that I can keep my step sizes fixed. This solver, however, blows up while the scipy LSODA (which adjusts for stiff equations) does not. Is there a name for this inherent stiffness or any literature you can recommend in order to study this further? $\endgroup$
    – Octavius
    Nov 27, 2022 at 14:44
  • $\begingroup$ I do not know if that has a special name somewhere. In a PDE $u_t=F(u,u_x,u_{xx},...)$ the partial derivatives on the right side are unbounded, smoothness-decreasing operators. The discretized versions approximate this property, the better the finer the grid. The inverse operators, resolvents, Green integral operators,... are inversely to the properties of the forward operator bounded, compact, smoothness-increasing. With implicit methods one can capture these properties for the integration step. (Which is the same story as for general ODE systems, only with a more geometric narrative.) $\endgroup$ Nov 27, 2022 at 14:55

2 Answers 2

4
$\begingroup$

If you use Julia's DifferentialEquations.jl, you can add adaptive=false to any method. solve(prob,DP5(),adaptive=false,dt=0.1) for RK45 and solve(prob,QNDF(),adaptive=false,dt=0.1) for a BDF method.

Now the reason why most solvers do not support this actually isn't too surprising. There's a lot of tricks that go into the behavior of Jacobian reuse which makes a lot of efficient adaptive methods for stiff ODEs fairly incompatible with just setting the dt to something fixed. The reason is that, if Newton methods diverge, the easiest thing to do is to simply decrease dt, snag a new Jacobian, and try again. Thus it is usually helpful to change the way the Newton method behaves if the adaptivity is removed. Sometimes setting always_new can help, i.e. change QNDF() to QNDF(nlsolve = NLNewton(check_div = false,always_new = true)), to force new Jacobians every time, like a naive university implementation.

$\endgroup$
1
  • $\begingroup$ Thank you for your answer! I will have a look into this Julia package you mentioned. $\endgroup$
    – Octavius
    Nov 27, 2022 at 14:47
4
$\begingroup$

I've answered this question a while ago on stackoverflow (I'll update with the link later).

Basically, you can get this with explicit schemes from solve_ivp (RK23, RK45, DOP853). I use this to impose a CFL condition for a compressible flow solver:

        dt = ... # Compute time step
        method.rtol=1e10
        method.max_step=dt
        method.first_step=dt
        method.h_abs=abs(dt)
        method.h    =abs(dt)

        method.step()
        t_new = method.t
        y_new = method.y

        dt_eff = t_new-t_old

        if not np.allclose(dt_eff, dt, rtol=1e-2, atol=1e-12):

            raise Exception('CFL was not maintained')

You may have missed some attributes of the method object (instance of the RK45 class for instance).

For implicit schemes (LSODA, BDF, Radau), this method does not work well, since the tolerance of the inner Newton loop is dependent on rtol... Thus convergence may be problematic. You would need to dig in the code to be able to prescribe an independent tolerance for the Newton solver.

$\endgroup$
2
  • $\begingroup$ Thank you for your answer. Could you elaborate shortly on the difference of method.h and method.h_abs? $\endgroup$
    – Octavius
    Nov 27, 2022 at 14:46
  • $\begingroup$ h may be negative if the final time is lower than the starting time in the tsoan argument. Otherwise, they have the same absolute value. I recommend reading through the step function of the RK methods of solve_ivp, everything will make sense ! $\endgroup$
    – Laurent90
    Nov 27, 2022 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.