1

I have an array of size (254, 80) which I need 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? Or perhaps if I am doing something wrong...

Here's the current code:

for i in range(array.shape[0]):
    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]) 

(The rest of the variables are all a similar size)

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?

Thanks, Rachael

1 Answer 1

0

Although noone answered I found a workaround which avoided the fsolve function and used interpolation instead. Luckily the initial guess is good enough that only a few y values are needed. If the initial guess knowledge is poor then this method is probably not appropriate. Do note this still has some issues but for my purposes it performs well...

ystar = np.empty((A,B))     # empty array for the solutions   
    num_ys = 20 #number of points to find where the solution is 
    y0_u = y0 #just so the calculated initial guess isn't overwritten
    for i in range(Thetas.shape[1]):
        ys = np.linspace(-.05,.2,num_ys)[:,None]*np.ones((num_ys,Thetas.shape[0])) + y0_u
        vals = (np.squeeze(eta) - np.squeeze(m)*ys*np.sqrt(g*np.tanh(ys**2*depth)) - np.squeeze(md)*np.sqrt(g*np.tanh(depth*np.sqrt(ys**4+2*(ys**2)*kB*np.cos(Thetas[:,i]+phi_bi)+kB**2)))*(( ys**4+2*(ys**2)*kB*np.cos(Thetas[:,i]+phi_bi)+kB**2 )**0.25))
        idxs_important = -1*(np.clip(np.vstack(((np.sign(vals[:-1]*vals[1:])-1),np.zeros((1,Thetas[:,i].size)))),-1,0) + np.clip(np.vstack((np.zeros((1,Thetas[:,i].size)),(np.sign(vals[:-1]*vals[1:]))-1)),-1,0))

        ys_chosen = idxs_important*ys
        ys_chosen[ys_chosen==0] = 10000 
        sorted_ys_idx = np.argsort(ys_chosen.T, axis = 1)
        sorted_ys = ((ys_chosen.T)[np.arange(np.shape(ys_chosen.T)[0])[:,np.newaxis],sorted_ys_idx]).T
        sorted_vals = (((vals*idxs_important).T)[np.arange(np.shape(vals.T)[0])[:,np.newaxis],sorted_ys_idx]).T
    #    interpolation bit 
        x_id = 0
        yposs = sorted_ys[:2,:]
        valposs = sorted_vals[:2,:]
        y = yposs[0,:] + (yposs[1,:] - yposs[0,:])*(x_id - valposs[0,:])/(valposs[1,:] - valposs[0,:])
        ystar[:,i] = np.squeeze(y)
        y0_u=ystar[:,i]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.