3

fsolve finds a solution of (a system of) nonlinear equations from a starting estimate. I can vectorize my function call to use fsolve on multiple starting points and potentially find multiple solutions, as explained here. In this question it is described how to solve multiple nonlinear equations with fsolve. However, I am having problems combining the two, i.e. solving multiple nonlinear equations from multiple starting values. I know that I can always loop through my starting values and use the second posts answer, but, having to do this for possibly more than a 100000 points, I'm really trying to find a more pythonic (and hopefully faster) solution.

I tried different ways, for example the following (and many others):

from scipy.optimize import fsolve
import numpy as np

def equations(x): # x+y^2-4, sin(x)+x*y-3
    ret = np.array([x[:,0]+x[:,1]**2-4, np.sin(x[:,0]) + x[:,0]*x[:,1] - 3]).T
    return ret

p1 = np.array([0,0]) # first initial value
p2 = np.array([1,1]) # second initial value
x0 = np.array([p1,p2])

print(x0[0,1])
print(equations(x0))
print(fsolve(equations, x0=x0))

The shapes and all work, but fsolve throws: 'IndexError: too many indices for array' I tried some different ways, but I can't work out any functioning code on this besides using a simple for loop. Any suggestions?

Banana
  • 1,149
  • 7
  • 24

2 Answers2

1

How about using joblib ? This is not a straight vectorization, but the different starting points will be executed in parallel.

from scipy.optimize import fsolve
import numpy as np
from joblib import Parallel, delayed

def equations(x): # x+y^2-4, sin(x)+x*y-3
    ret = np.array([x[0]+x[1]**2-4, np.sin(x[0]) + x[0]*x[1] - 3]).T
    return ret

p1 = np.array([0,0]) # first initial value
p2 = np.array([1,1]) # second initial value
x0 = [p1, p2]

sol = Parallel(n_jobs=2)(delayed(fsolve)(equations, x0=p) for p in x0)
print(sol)

The parameter n_jobs controls how many concurrent jobs are running.

fakufaku
  • 96
  • 6
0
def eq(x):
     return x[0] + x[1]**2 - 4 , np.sin(x[0]) + x[0]*x[1] - 3

fsolve(eq, [0, 1])

output: array([ 1.23639399,  1.6624097 ])

In this case, I would suggest using a brute force method:

x0 = [[i, j] for i, j in zip(range(10), range(10))]

for xnot in x0:
    print(fsolve(eq, xnot))

[  1.33088471e-07   2.09094320e+00]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
[ 1.23639399  1.6624097 ]
Osman Mamun
  • 2,864
  • 1
  • 16
  • 22