Faster method to numerously solve for multiple parameters governed by multiple nonlinear equations with multiple variable arguments?

I am looking for a tweak that I could apply to the code shown later or an alternate method that would lead to faster run time. In general, I want to find the values of two parameters that are governed by a set of multiple nonlinear equations. These equations use multiple variable arguments.

An Example of the problem at hand is presented for only a small array (4 values) as an example:

I want to find the values of q and h that are governed by 2 non-linear equations
I will need to apply this procedure to larger a and b arrays each of about 50k values

Equations are:

q = k/2L * (h^2 - b^2)
q = -2 * pi * k * (h-a) / ln((h-d)/r)

arguments/variables:

k = 144.0
L = 550.0
d = 140.9
r = 0.5
a = [190.0, 185.0, 160.0, 150.0]
b = [70.0, 70.0, 30.0, 10.0]

I followed suggestions in similar posts at stackoverflow, mostly This post, and modified my original code to avoid using a for loop. Though, after doing that, I tried applying the code to a big array of variables, but it took quite some time.

I should also mention that I am only interested in the solution when h is between a certain limit (when h is greater than d+r and less than 300), but I didn't know how to apply this in fsolve when I am solving for multiple parameters and using multiple non-linear functions.

The final code that I am currently using is shown in this post, but I am looking for something that is faster and more efficient.

import numpy as np

from scipy import optimize

def functions(X, Y):

    q,h = np.split(X,2)

    a,b = np.split(Y[:-4],2)

    k,L = Y[-4],Y[-3]

    d,r = Y[-2],Y[-1]

    f1 = q - k/(2.0*L)*(h**2 - b**2) 

    f2 = q + 2.0*np.pi*k*(h - a)/ np.log((h - d)/r)

    return np.append(f1, f2)

k = 144.0

L = 550.0

d = 140.9

r = 0.5

a_vals = np.array([190.0, 185.0, 160.0, 150.0])

b_vals = np.array([70.0, 70.0, 30.0, 10.0])

q0 = np.array([4059.5, 3814.7, 3212.5, 2912.7]) #initial guesses for q

h0 = np.array([189.5, 184.5, 159.5, 149.5]) #initial guesses for h

my_args = np.concatenate((a_vals, b_vals, k, L, d, r), axis=None)

initial_guesses = np.concatenate((q0, h0), axis=None)

qs,hs = np.split(optimize.fsolve(func = functions, x0 = initial_guesses, args= my_args), 2)

Solution

The simplest solution in your case is to use the Numba's JIT since many (short) calls to the same function are performed. Here is an example:

import numba as nb

@nb.njit('float64[::1](float64[::1],float64[::1])')
def functions(X, Y):
    q,h = np.split(X,2)
    a,b = np.split(Y[:-4],2)
    k,L = Y[-4],Y[-3]
    d,r = Y[-2],Y[-1]
    f1 = q - k/(2.0*L)*(h**2 - b**2) 
    f2 = q + 2.0*np.pi*k*(h - a)/ np.log((h - d)/r)
    return np.append(f1, f2)

# ... (the rest is the same)

This is 12 times faster on my machine.