2

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.

J.A
  • 171
  • 3
  • 2
    What have you already done? It's just a regular ODE, for which there are a large number of step length choice methods. – Wolfgang Bangerth Apr 01 '21 at 13:41
  • Going off of @WolfgangBangerth comment my first reaction would be to look at the eigenvalues of $a_{ij}$ and make sure they are in the stability region of forward Euler. – Kyle Mandli Apr 01 '21 at 20:00
  • @KyleMandli This would be appropriate if the equation was linear, but it is quadratic in the right hand side. One needs to look at the eigenvalues of $a_{ij}y_i$. – Wolfgang Bangerth Apr 01 '21 at 22:50
  • @WolfgangBangerth I added what I tried – J.A Apr 02 '21 at 08:17
  • @KyleMandli all the eigenvalues must be negative right ? But does this provide a time-step ? – J.A Apr 02 '21 at 08:19
  • 2
    The eigenvalues of the discretized system must be stable, so for forward euler, you must achieve $|1 + \lambda \delta t| < 1$ – Thijs Steel Apr 02 '21 at 08:28
  • 4
    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. – Laurent90 Apr 02 '21 at 14:07
  • 1
    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. – Kyle Mandli Apr 07 '21 at 00:24
  • @KyleMandli I tried with odeint. But it is quite slow... Are there other adaptive time-steppers ? – J.A Apr 08 '21 at 09:44
  • 1
    @J.A Use the integrators from scipy.integrate.solve_ivp. A good explicit integrator is 'RK45', equivalent to ode45in Matlab. Play with the integration tolerances rtol and atol(set them to 1e-3, 1e-6 or 1e-9 for instance) and see how well (or not) your problem is solved. If the solver is slow, I recommend trying the implicit integrators 'BDF' (or LSODA which is nearly identical) and 'Radau'. Ideally, post a standalone code snippet that allows us to easily run your test case and play with it ;) – Laurent90 Apr 08 '21 at 10:25
  • 2
    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. – Lutz Lehmann Apr 08 '21 at 11:48

0 Answers0