调试，solve_ivp与toms748一起使用时

Question

我在执行以下代码时遇到错误

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

得到的错误如下

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 是一个求根算法，它接受一个可调用函数和两个标量，这两个标量绑定了要在其间搜索根的值。因此，E(t,r) 就是 E(t,r)=-r，上面实现的微分方程是 dr/dt=-r，初始条件为 r(0)=1。解决方案只是 r(t)=exp(-t).

现在更让我困惑的是，当我从 E(t,r) 中删除减号时，即让

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

然后代码运行正常，没有抛出任何错误。

以上都是玩具功能。实际实现要复杂得多。上面是我能得到的最简单的代码，它给出了同样的错误。

Answer 1

toms748 期望使用标量，而您给它一个长度为 1 的向量。如果我将 E 更改为

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

ODE 求解器运行完成。

为了帮助调试，我将您的代码更改为：

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

这给了我：

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

我可以单独使用 E 重现错误，使用 r.

的最后一个全精度值

如果toms748能给出更好的错误信息就更好了，我不知道为什么有些调用失败有些成功——大概代码中的其他分支使用了不同的方式来索引。

调试，solve_ivp与toms748一起使用时

Debugging, solve_ivp when used with toms748

python

scipy

ode