有没有一种方法可以将fsolve矢量化？

Question

我想申请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 

解决它。

Answer 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。

Answer 2

可以定义AA功能，以尽量减少（这应该是你方原有的功能），然后用一个简单的极小（你最好也定义函数的导数）：

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)