I have a system of non linear differential equations in the form : $$\frac{dy_i}{dt}=\sum_j a_{ij} y_i y_j $$.

I first tried to solve it with Python suing scipy.integrate.odeint but it is very slow.

I'd like to solve it with an Euler method : $$y_i(t+\delta t)=y_i(t)+\frac{dy_i}{dt}\delta t $$.

How should I chose the value of $\delta t$ to ensure a good convergence ?

What I tried is the following.

Let's define the error $\epsilon$ : $$\epsilon=\sqrt{\sum_i(\delta t \times dy_i/dt)^2} $$

Starting from a low enough initial guess for $\delta t$, I looked wether the error increase or descrease.

  • If the error decreases : $\epsilon(t)>\epsilon(t+\delta t)$, I change $\delta t\rightarrow 2\delta t$
  • If the error increases : $\epsilon(t)<\epsilon(t+\delta t)$, I change $\delta t\rightarrow \delta t/2$

But the calculation is very slow and can even diverge.

Another possibility I tried :

  • If the error decreases : $\epsilon(t)>\epsilon(t+\delta t)$, I change $\delta t\rightarrow 2^n\delta t$ with $n$ so that $\epsilon(t+2^{n}\delta t)<\epsilon(t)<\epsilon(t+2^{n+1}\delta t)$
  • If the error increases : $\epsilon(t)<\epsilon(t+\delta t)$, I change $\delta t\rightarrow \delta t/2^n$ until $\epsilon(t+\delta t/2^{n})<\epsilon(t)<\epsilon(t+\delta t/2^{n-1})$

but this scheme is not very efficient too.

  • 2
    $\begingroup$ What have you already done? It's just a regular ODE, for which there are a large number of step length choice methods. $\endgroup$ Apr 1, 2021 at 13:41
  • 2
    $\begingroup$ The eigenvalues of the discretized system must be stable, so for forward euler, you must achieve $|1 + \lambda \delta t| < 1$ $\endgroup$ Apr 2, 2021 at 8:28
  • 4
    $\begingroup$ Note that this stability criterion does not guarantee that the solution obtained is accurate. For that, you may want to look into adaptive time stepping (embedded Runge-Kutta method among others). For forward Euler, you can create, using the solution values/derivatives at the previous time step, a 2nd order approximation of your solution via polynomial interpolation (as for multistep methods), from which you can form an estimate of the local error to determine a suitable time step. More details can be found in any book on numerical methods for ODEs. $\endgroup$
    – Laurent90
    Apr 2, 2021 at 14:07
  • 1
    $\begingroup$ Have you tried one of the basic adaptive time steppers from something from SciPy or even PETSc? They may be able to handle the non-linearity you are seeing. $\endgroup$ Apr 7, 2021 at 0:24
  • 2
    $\begingroup$ If odeint is slow, and similarly solve_ivp with method="Radau", then the system is stiff. Which is not surprising for a quadratic right side. It would, on the other hand, require some structure in the coefficients to prevent the solution from blowing up. You will not get a useful solution faster with the Euler method. You would need step sizes that are so small that the effort becomes much higher. // Due to the algebraic construction of the right side, you might be able to employ operator splitting methods, like the exponential variants of Runge-Kutta. $\endgroup$ Apr 8, 2021 at 11:48


Your Answer

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

Browse other questions tagged or ask your own question.