Есть ли способ векторизовать fsolve?

Я пытаюсь применить 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
person Adobe    schedule 31.05.2012    source источник

Ответы (2)


arrow_upward
3
arrow_downward

fsum предназначен для скаляров Python, поэтому вам следует обратиться к numpy для векторизации. Ваш метод, вероятно, терпит неудачу, потому что вы пытаетесь суммировать список из пяти массивов numpy, а не пяти чисел или одного массива numpy.

Сначала я бы пересчитал cn с помощью numpy:

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

or

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

fsolve, кажется, подбирает первое решение при запуске с a0 = 0, а второе - для a0 = 3.

person marshall.ward    schedule 31.05.2012
comment
Бит это дает вам один a, в то время как их должно быть 400 - для каждого b. Смотрите редактирование исходного поста. - person Adobe; 31.05.2012
comment
Вы правы, мне нужно было использовать вектор в качестве начального предположения. Я обновил это предложение. - person marshall.ward; 31.05.2012

arrow_upward
1
arrow_downward

Вы можете определить функцию для минимизации (которая должна быть квадратом вашей исходной функции), а затем использовать простой минимизатор (вам также лучше определить производную функции):

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)
person sega_sai    schedule 31.05.2012