CODEWITHSUNDEEP
pythonscipydifferential-equations
PingPong
PingPong

Reputation: 41

Solving ODE with solve_ivp gets incredibly slow or freezes completely

I use solve_ivp to solve an ODE:

def test_ode(t, y):
    dydt = C - y + (y ** 8 / (1 + y ** 8))
    return dydt

steady_state = []
for C in np.linspace(0, 1, 1001):
    sol = solve_ivp(test_ode, [0, 1e06], [0], method='BDF')
    steady_state.append(sol.y[0][-1])

This gives me a RuntimeWarning:

ETA:  --:--:--/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/bdf.py:418: 
RuntimeWarning: divide by zero encountered in power
  factors = error_norms ** (-1 / np.arange(order, order + 3))

But even worse, the run basically freezes (or at least gets incredibly slow). Replacing the initial value [0] by [1e-08] does not solve the problem. How can I fix this?

Upvotes: 1

Views: 3588

Answers (1)

nicholaswogan
nicholaswogan

Reputation: 679

You can use NumbaLSODA: https://github.com/Nicholaswogan/NumbaLSODA . Its like solve_ivp, but all the code can be compiled. So it is very speedy:

from NumbaLSODA import lsoda_sig, lsoda
import numpy as np
import numba as nb
import time

@nb.cfunc(lsoda_sig)
def test_ode(t, y_, dydt, p):
    y = y_[0]
    C = p[0]
    dydt[0] = C - y + (y ** 8 / (1 + y ** 8))
funcptr = test_ode.address

@nb.njit()
def main():
    steady_state = np.empty((1001,),np.float64)
    CC = np.linspace(0, 1, 1001)
    y0 = np.array([0.0])
    for i in range(len(CC)):
        data = np.array([CC[i]],np.float64)
        t_eval = np.array([0.0, 1.0e6])
        sol, success = lsoda(funcptr, y0, t_eval, data = data)
        steady_state[i] = sol[-1,0]
    return steady_state

main()
start = time.time()
steady_state = main()
end = time.time()
print(end-start)

result is 0.013 seconds

Upvotes: 1

Related Questions