有沒有一種方法可以將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)