# solve_ivp doesn't work with toms748

I have 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])


When I run it I get the following error

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-
_, 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):


then no errors are returned

All of the above are toy functions. The actual functions are much more complicated and the above is just the bare bones implementation that reproduces the same errors. Any help is greatly appreciated, thanks in advance.

I'm not certain of the exact error, but the issue seems to stem from tom748 when you pass in an array/list r to E rather than a scalar, though I'm not sure why it only seems to happen when solveivp uses -tom748.

solveivp requires you to enter your initial value as an array-like, but you can rewrite E to unpack the element from the list:

def E(t,r):
r0=r[0]
res= -toms748(f(r0),r0-1,r0+1)
return res


With this modification, the script ran just fine for me.

• First of all thanks a bunch. So it seems like even if some r values are passed as arrays in toms748 then no error is triggered. For instance if we consider the same f(r) as defined above then toms748(f(2), np.array([-10]), np.array([10])) gives an error but not if the first bound is replaced by np.array([-1]). Therefore my guess is that when the positive value for toms748 is used these ''trigger'' values are not reached in the course of integration. Nov 20, 2021 at 18:00