Debugging, solve_ivp when used with toms748

Question

I am facing an error when implementing the following code

from scipy.optimize import toms748
from scipy.integrate import solve_ivp

def f(r):
    return lambda x: x-r

def E(t,r):
    return -toms748(f(r),r-1,r+1)

sol=solve_ivp(E,(0,10),[1])

The error obtained is as follows

Traceback (most recent call last):
File "C:\Users\User\Documents\Project codes\Density.py", line 10, in <module>
sol=solve_ivp(E,(0,10),[1])
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\integrate\_ivp\ivp.py", line 576, in solve_ivp
message = solver.step()
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\integrate\_ivp\base.py", line 181, in step
success, message = self._step_impl()
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\integrate\_ivp\rk.py", line 144, in _step_impl
y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\integrate\_ivp\rk.py", line 64, in rk_step
K[s] = fun(t + c * h, y + dy)
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\integrate\_ivp\base.py", line 138, in fun
return self.fun_single(t, y)
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\integrate\_ivp\base.py", line 20, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
File "C:\Users\User\Documents\Project codes\Density.py", line 8, in E
return -toms748(f(r),r-1,r+1)
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\optimize\zeros.py", line 1361, in toms748
result = solver.solve(f, a, b, args=args, k=k, xtol=xtol, rtol=rtol,
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\optimize\zeros.py", line 1225, in solve
status, xn = self.iterate()
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\optimize\zeros.py", line 1144, in iterate
c = _newton_quadratic(self.ab, self.fab, d, fd, nsteps)
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\optimize\zeros.py", line 1004, in _newton_quadratic
_, B, A = _compute_divided_differences([a, b, d], [fa, fb, fd],
File "C:\Users\User\AppData\Local\Programs\Python\Python39\lib\site- 
packages\scipy\optimize\zeros.py", line 959, in _compute_divided_differences
row = np.diff(row)[:] / denom
ValueError: operands could not be broadcast together with shapes (3,0) (2,1)

toms748 is a root finding alogrithm that takes in a callable function and two scalars that bound the values between which the root is searched for. Thus E(t,r) is just E(t,r)=-r and the differential equation implemented above is dr/dt=-r with initial condition r(0)=1. The solution is just r(t)=exp(-t).

Now the thing that perplexes me even more is, when I remove the minus sign from E(t,r) i.e let

def E(t,r):
    return toms748(f(r),r-1,r+1)

then the code runs just fine and no errors are thrown up.

The above are all toy functions. The actual implementation is much more complicated. The above is the simplest code I could obtain that gives the same error.

I've looked into the solver code in connection with a different question, see https://stackoverflow.com/questions/59389223/modelling-a-quantum-harmonic-oscillator. I do not know if the behavior has changed in-between. My impression was that the `toms748` is legacy code, recommended is to use `root_scalar` and, if necessary, set the method parameter for this procedure. — Lutz Lehmann, Nov 20 '21 at 16:49

score 1 · Accepted Answer · answered Nov 20 '21 at 15:18

toms748 is expecting to work with scalars, and you're giving it a length-1 vector. If I change E to

def E(t,r):
    print(f'{t=!r}, {r=!r}, {r[0]:.17g}')
    return -toms748(f(r[0]),r[0]-1,r[0]+1)

the ODE solver runs to completion.

To help debug this, I changed your code to:

from scipy.optimize import toms748
from scipy.integrate import solve_ivp

def f(r):
    def myf(x):
        print(f'{x=!r}, {x-r=!r}')
        
        return x-r
    return myf

def E(t,r):
    print(f'{t=!r}, {r=!r}, {r[0]:.17g}')
    return -toms748(f(r),r-1,r+1)

sol=solve_ivp(E,(0,10),[1])

That gave me:

t=0.0, r=array([1.]), 1
x=array([0.]), x-r=array([-1.])
x=array([2.]), x-r=array([1.])
x=array([1.]), x-r=array([0.])
t=0.01, r=array([0.99]), 0.98999999999999999
x=array([-0.01]), x-r=array([-1.])
x=array([1.99]), x-r=array([1.])
x=array([0.99]), x-r=array([0.])
t=0.020003998400959323, r=array([0.979996]), 0.97999600159904066
x=array([-0.020004]), x-r=array([-1.])
x=array([1.979996]), x-r=array([1.])
x=array([0.979996]), x-r=array([1.11022302e-16])

and I could reproduce the error with E alone using that last full-precision value for r.

It would have been ideal if toms748 would give a better error message, and I don't know why some calls fail and some succeed - presumably other branches in the code use a different means to index.

Thanks a ton. Really appreciate writing out the debugging/thought process. So it seems like even when `toms748` is given size 1 arrays the code runs for a while and then some value of `r` trips it up. Just for curiosity's sake do you have any idea why only certain values would affect `toms748` when arrays are given as inputs? — user7896, Nov 21 '21 at 05:11
No - as I said, my guess is different branches in `toms78` use different indexing methods. This sort of value-dependent error can be a pain to debug. — Rory Yorke, Nov 21 '21 at 05:57

Debugging, solve_ivp when used with toms748

1 Answers1