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 '21 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 '21 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 '21 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 '21 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 '21 at 11:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.