Find time step for Euler method in numerical solving of a system of non linear differential equations

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.

• What have you already done? It's just a regular ODE, for which there are a large number of step length choice methods. Apr 1 '21 at 13:41
• The eigenvalues of the discretized system must be stable, so for forward euler, you must achieve $|1 + \lambda \delta t| < 1$ Apr 2 '21 at 8:28
• 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. Apr 2 '21 at 14:07
• 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. Apr 7 '21 at 0:24
• 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. Apr 8 '21 at 11:48