0

Let us assume I have an ODE with x'(t) = f(x) with the respective solution x(t) = ϕ(x(0),t) of a initial condition x(0). Now I intend to calculate numerically the equilibria as a function of their initial condition: eq(x0) := ϕ(x0, ∞). The ODEs are such that these equilibria exist unambiguously for all initial conditions (including eq = ∞).

My poor man's approach would be to integrate the ODE up to a late time and fetch that value (for brevity I do not show the plotting):

import numpy as np
from scipy.integrate import odeint

# ODE
def func(X,t):
    return [ X[2]**2 * (X[0] - X[1]),
             X[2]**3 * (X[0] + 3 * X[1]),
             -X[2]**2] 

# Forming a grid
n = 15
x0 = x1 = np.linspace(0,1,n)
x0_,x1_ = np.meshgrid(x0,x1)
eq = np.zeros([n,n,3])

t  = np.linspace(0,100,1000)  
x2 = 1

for i in range(n):
    for j in range(n):
        X = odeint(func,[x0_[j,i],x1_[j,i],x2], t)
        eq[j,i,:] = X[-1,:]

Naive example above:

enter image description here

The problem with that approach is that you can never be sure if it converged. I know that you can just find the roots of f(x), but this would not yield the equilibria as a function of their initial conditions (You could trace them back, but since this function is not injective, you will not find values for all initial values). I somehow need a ODE solver which integrates until an equilibria is reached (or stops integrating if it goes beyond a limit). Do you have any ideas?

varantir
  • 6,624
  • 6
  • 36
  • 57
  • I consider this question off-topic here as it is asking for algorithmic approaches and has nothing to do with programming. It is probably best suited for [scicomp.se]. – Wrzlprmft Sep 08 '17 at 12:25
  • Acutally I am searching for an alternative to scipy.odeint since in this scenario it is not possible to introduce a stop-condition for the integration. I would consider this to be a programming-related question – varantir Sep 08 '17 at 13:40
  • I thought your problem was to find some criterion for stopping the integration in the first place. If you already have that the criterion, you can strongly boil down your question. – Wrzlprmft Sep 08 '17 at 13:59
  • 1
    From question [Can I integrate with scipy's odeint until a local max is found?](https://stackoverflow.com/q/30597976/5358968), it looks like you can call `integrate` to a time, test, then integrate further. From [Does scipy.integrate.ode.set_solout work?](https://stackoverflow.com/q/26738676/5358968), it looks like you can provide a stopping condition with `set_solout`. I didn't get much value from looking for zeroes since I got muddled with the `x_2` term in `f_1`. – Steve Sep 11 '17 at 09:23

0 Answers0