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.
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.