使用循环求根

Question

我在函数中定义了一个方程

def fun(x, y, z, v, b):
       Y = (z*(np.sign(x) * (np.abs(x))**(y-1))) - (v*np.sign(b) * (np.abs(b))**(v-1))/(1-b**v)
       return Y.flatten()

我想求解 x 的值，给定 Z0、SS（第 1 年：Z0=1.2，SS=2，...）的值以及 alpha 和 kappa，我正在为其创建一个网格。

Z0 = [1.2, 5, 3, 2.5, 4.2]
SS = [2, 3, 2.2, 3.5, 5]

ngrid = 10
kv = np.linspace(0.05, 2, ngrid)
av = np.linspace(1.5, 4, ngrid)

q0 = []                    
for z in range(len(Z0)):
    zz = Z0[z]
    ss = SS[z]
    for i in range(ngrid):
        for j in range(ngrid):
            kappa = kv[i]
            alpha = av[j]
            res0 = root(lambda x: fun(x, alpha, zz, kappa, ss), x0=np.ones(range(ngrid)))
            q0 = res0.x
            print(q0)

其中 y = α； v=kappa，z=Z0； b = S.

我得到了所有 [], [], ....

不确定发生了什么。感谢您的帮助

Answer 1

在您尝试使用 res0.x 之前，请检查 res0.success。在这种情况下，您会发现它在每种情况下都是 False。当 res0.success 为 False 时，请查看 res0.message 了解有关 root 失败原因的信息。

在开发和调试期间，在将 root 嵌入三个嵌套循环之前，您还可以考虑让求解器仅处理一组参数值。例如，这里是 ipython 会话中的几行（变量在前面的行中定义，未显示）：

In [37]: res0 = root(lambda x: fun(x, av[0], Z0[0], kv[0], SS[0]), x0=np.ones(range(ngrid)))

In [38]: res0.success
Out[38]: False

In [39]: res0.message
Out[39]: 'Improper input parameters were entered.'

消息提示输入参数有问题。你这样调用 root：

        res0 = root(lambda x: fun(x, alpha, zz, kappa, ss), x0=np.ones(range(ngrid)))

仔细看那一行就会发现问题所在：最初的猜测是 np.ones(range(ngrid)):

In [41]: np.ones(range(ngrid))
Out[41]: array([], shape=(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), dtype=float64)

这不是你想要的！ range 的使用看起来像是一个简单的错字（或“thinko”）。初步猜测应该是

x0=np.ones(ngrid)

在ipython中，我们得到：

In [50]: res0 = root(lambda x: fun(x, av[0], Z0[0], kv[0], SS[0]), x0=np.ones(ngrid))

In [51]: res0.success
Out[51]: True

In [52]: res0.x
Out[52]: 
array([-0.37405428, -0.37405428, -0.37405428, -0.37405428, -0.37405428,
       -0.37405428, -0.37405428, -0.37405428, -0.37405428, -0.37405428])

所有 return 值都相同（其他参数值也会出现这种情况），这表明您正在求解标量方程。仔细查看 fun 表明您仅在逐元素运算中使用 x，因此您实际上只是在求解一个标量方程。在这种情况下，您可以使用 x0=1:

In [65]: res0 = root(lambda x: fun(x, av[0], Z0[0], kv[0], SS[0]), x0=1)

In [66]: res0.success
Out[66]: True

In [67]: res0.x
Out[67]: array([-0.37405428])

您也可以考虑使用 root_scalar 而不是 root。

使用循环求根

Root finding using a loop

root

scipy