Is there a way I can vectorize fsolve?

Question

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)

It is not allowed by default:

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

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

Edit:

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

solves it.

marshall.ward · Accepted Answer · 2012-05-31T12:05:04.297

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.

First I would recalculate cn using numpy:

import numpy as np

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

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)

Redefining your function in terms of K should now work:

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

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)

or you can use a0 = 3:

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

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

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

or

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

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

Bit it gives You a single `a` while it should be 400 of them - for each `b`. See the edit of the original post. — Adobe, May 31 '12 at 11:35
You are right, I needed to use a vector as the initial guess. I've updated this suggestion. — marshall.ward, May 31 '12 at 12:02

sega_sai · Answer 2 · 2012-05-31T12:41:34.623

You can define a a function to minimize (which should be a square of you original function) and then use a simple minimizer (you better also define the derivative of the function):

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)

Is there a way I can vectorize fsolve?

2 Answers2

Linked