I have an array of size (254, 80) which I am trying to use Scipy's fsolve on. I have found that the speed of using fsolve on a vector is quicker than it is in a for loop but only for vectors upto about 100 values long. After this, the speed quickly drops off and becomes very slow, sometimes completely stopping.
I'm currently looping through one dimension of the array and using a vectorised fsolve on the smaller dimension but it's still taking longer than I would expect/like.
Does anyone have a good work around for this or a know of a similar function which will be happy handling a vector of a larger size? I have tried Scipy's Newton Krylov function which, upon testing, is really fast for simple functions but it's slower than fsolve on my function. - Will this be faster if I use the analytical jacobian? If so, how do I implement it?
I wondered if it is a problem with my initial solutions but the difference is in the range of ~0.001 - 0.8 which doesn't seem too bad. Although perhaps it is this?
Here's the current code:
for i in range(array.shape): f = lambda y: a[i] - m[i]*y - md[i]*(( y**4 + 2*(y**2)*np.cos(Thetas[i,:]) )**0.25) ystar[i,:] = fsolve(f, y0[i])
(Variables a, m and md are all the same size i.e. (254,80))
Digging in to this further, I have found that a function such as
f = lambda y: y*np.tanh(y) - a0/(m**2)
is faster to solve than
f = lambda y: (m**2)y*np.tanh(y) - a0
where m and a0 are large 2D np arrays.
Can anyone explain why this is?