有没有办法可以矢量化 fsolve? [英] Is there a way I can vectorize fsolve?

问题描述

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

I'm trying apply fsolve to an array:

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 应用于数组?

is there a way I can apply fsolve to an array?

编辑:

#!/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

解决它.

推荐答案

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

fsum is for python scalars, so you should look to numpy for vectorisation. Your method is probably failing because you're trying to sum a list of five numpy arrays, rather than five numbers or a single numpy array.

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

First I would recalculate cn using numpy:

import numpy as np

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

接下来我将分别计算之前的 fsum 结果，因为它是一个常数向量.这是一种方式，虽然可能有更有效的方式:

Next I would compute the previous fsum result separately, since it's a constant vector. This is one way, although there may be more efficient ways:

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

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

Redefining your function in terms of K should now work:

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

<小时>

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

---

To use fsolve to find the solution, provide it with an appropriate initial vector to iterate against. This uses the zero vector:

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

或者你可以使用 a0 = 3:

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

这个函数是可逆的，所以你可以检查 f(a) = 0 对两个确切的解决方案:

This function is invertible, so you can check f(a) = 0 against the two exact solutions:

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

或

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

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

fsolve seems to be picking up the first solution when starting from a0 = 0, and the second one for a0 = 3.

这篇关于有没有办法可以矢量化 fsolve?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！