有没有办法可以矢量化 fsolve？答案

【问题标题】：Is there a way I can vectorize fsolve?有没有办法可以矢量化 fsolve？
【发布时间】：2012-05-31 20:41:34
【问题描述】：

我正在尝试将 fsolve 应用于数组：

from __future__ import division
from math import fsum
from numpy import *
from scipy.optimize import fsolve
from scipy.constants import pi

nu = 0.05
cn = [0]
cn.extend([pi*n - pi/4 for n in range(1, 5 +1)])
b = linspace(200, 600, 400)
a = fsolve(lambda a: 1/b + fsum([1/(b+cn[i]) for i in range(1, 5 +1)]) - nu*(sqrt(a) - 1)**2, 3)

默认是不允许的：

TypeError: only length-1 arrays can be converted to Python scalars

有没有办法可以将 fsolve 应用于数组？

编辑：

#!/usr/bin/env python

from __future__ import division
import numpy as np
from scipy.optimize import fsolve

nu = np.float64(0.05)

cn = np.pi * np.arange(6) - np.pi / 4.
cn[0] = 0.

b = np.linspace(200, 600, 400)

cn.shape = (6,1)
cn_grid = np.repeat(cn, b.size, axis=1)
K = np.sum(1/(b + cn_grid), axis=0)

f = lambda a: K - nu*(np.sqrt(a) - 1)**2
a0 = 3. * np.ones(K.shape)
a = fsolve(f, a0)

print a

解决了。

【问题讨论】：

标签： numpy scipy

【解决方案1】：

fsum 用于 python 标量，因此您应该使用 numpy 进行矢量化。您的方法可能失败了，因为您尝试对五个 numpy 数组的列表求和，而不是五个数字或单个 numpy 数组。

首先我会使用 numpy 重新计算 cn：

import numpy as np

cn = np.pi * np.arange(6) - np.pi / 4.
cn[0] = 0.

接下来我将分别计算前一个 fsum 结果，因为它是一个常数向量。这是一种方法，尽管可能有更有效的方法：

cn.shape = (6,1)
cn_grid = np.repeat(cn, b.size, axis=1)
K = np.sum(1/(b + cn_grid), axis=0)

根据K 重新定义您的函数现在应该可以工作了：

f = lambda a: K - nu*(np.sqrt(a) - 1)**2

要使用fsolve 找到解决方案，请为其提供适当的初始向量以进行迭代。这使用了零向量：

a0 = np.zeros(K.shape)
a = fsolve(f, a0)

或者你可以使用a0 = 3:

a0 = 3. * np.ones(K.shape)
a = fsolve(f, a0)

这个函数是可逆的，所以你可以检查f(a) = 0 和两个精确的解：

a = (1 - np.sqrt(K/nu))**2

或

a = (1 + np.sqrt(K/nu))**2

fsolve 从a0 = 0 开始时似乎选择了第一个解决方案，而a0 = 3 则选择了第二个解决方案。

【讨论】：

位它给你一个a，而它应该是400个——每个b。请参阅原始帖子的编辑。
你说得对，我需要使用向量作为初始猜测。我已更新此建议。

【解决方案2】：

你可以定义一个函数来最小化（它应该是你原始函数的平方），然后使用一个简单的最小化器（你最好也定义函数的导数）：

funcOrig = lambda a: (K - nu*(np.sqrt(a) - 1)**2)
func2 = lambda a: funcOrig(a)**2
dotfunc2 = lambda a: 2*funcOrig(a)* (-nu * 2 * ( np.sqrt(a)-1) * 1./2./np.sqrt(a))
ret = scipy.optimize.fmin_l_bfgs_b(func2, np.ones(400)+1, fprime=dotfunc2, pgtol=1e-20)

【讨论】：