Fastest way to solve an array or list of functions with fsolve

user_3173758 0
Viewed 267

I have the working function below. I have a function from which I calculate the first and second derivative. I then need to find the value of theta where the first derivative is zero and second one is negative. I have to compute this for a large number of points. The number of points is equal to the length of K1 and K2. Using sympy I calculate the first and second derivative. I currently iterate over all the derivatives and solve the equations for each. Is there a faster way to do this, once the lengths of K1 and K2 increase > 1000 , this takes much to long for my application.

import numpy as np
import sympy as sp
from scipy.optimize import fsolve
from sympy.utilities.lambdify import lambdify

def get_cpd(K1, K2):
       '''
        *args:
            K1, 1D numpy array: mode I stress intensity factors
            K2, 1D numpy array: mode II stress intensity factor
        *Note:
            K1 and K2 should have the same length
       '''
        # Define symbols
        r, theta = sp.symbols("r theta")

        # Shear stress intensity
        sif_shear = 1/2*sp.cos(theta/2)*(K1*sp.sin(theta)+K2*(3*sp.cos(theta)-1))

        # Determine the first and second derivative w.r.t. theta
        first_derivative = sp.diff(sif_shear, theta)
        second_derivative = sp.diff(first_derivative, theta)

        cpd_lst = []
        for first, second in zip(first_derivative, second_derivative):
            # Lambdify function such that it can evaluate an array of points
            func1 = sp.lambdify(theta, first, "numpy")
            func2 = sp.lambdify(theta, second, "numpy")

            # initialize array from -π/2 to π/2, this is used for the initial guesses of the solver
            x = np.linspace(-np.pi/2, np.pi/2, num=50)

            # Solve the first derivative for all initial guesses to find possible propagation angles
            y1 = fsolve(func1, x)

            # Evaluate the second derivative in the roots of the first derivative
            y2 = func2(y1)
                   
            # Get roots of first derivative between -π/2 to π/2 
            # and where second derivative is negative
            y1 = np.round(y1, 4)
            y1 = y1[(y1 > -np.pi/2) & (y1 < np.pi/2) & (y2 < 0)]

            # get unique roots
            cpd = np.unique(y1)

            cpd_lst.append(cpd)

        return cpd_lst

Input example:

K1 = np.random.rand(10000,)
K2 = np.random.rand(10000,)
get_cpd(K1, K2)
sympy scipy numpy python
1 Answers

The best thing is to try and process the equation symbolically as much as possible in terms of symbolic parameters. It is possible to get an analytic solution for e.g. first_derivative but you need to transform it a bit. Here I'll rewrite sin/cos as exp and then use the substitution exp(I*theta/2) = sqrt(z) to get a cubic polynomial for z:

In [150]: K1, K2 = symbols('K1, K2', real=True)                                                                                                                                                             

In [151]: theta = Symbol('theta', real=True)                                                                                                                                                                

In [152]: sif_shear = S.Half*sp.cos(theta/2)*(K1*sin(theta)+K2*(3*cos(theta)-1))                                                                                                                            

In [153]: eq = diff(sif_shear, theta)                                                                                                                                                                       

In [154]: eq                                                                                                                                                                                                
Out[154]: 
  ⎛K₁⋅sin(θ)   K₂⋅(3⋅cos(θ) - 1)⎞    ⎛θ⎞                                   
  ⎜───────── + ─────────────────⎟⋅sin⎜─⎟                                   
  ⎝    2               2        ⎠    ⎝2⎠   ⎛K₁⋅cos(θ)   3⋅K₂⋅sin(θ)⎞    ⎛θ⎞
- ────────────────────────────────────── + ⎜───────── - ───────────⎟⋅cos⎜─⎟
                    2                      ⎝    2            2     ⎠    ⎝2⎠

In [155]: eqz = fraction(cancel(eq.rewrite(exp).subs(exp(I*theta/2), sqrt(z))))[0].collect(z)                                                                                                               

In [156]: eqz                                                                                                                                                                                               
Out[156]: 
                 3                    2                            
3⋅K₁ - 9⋅ⅈ⋅K₂ + z ⋅(3⋅K₁ + 9⋅ⅈ⋅K₂) + z ⋅(K₁ + ⅈ⋅K₂) + z⋅(K₁ - ⅈ⋅K₂)

Now sympy can solve this (roots(eqz, z)) but the general formula for a cubic is quite complicated so that might not be the best approach. Given particular float values for K1 and K2 though sympy can easily get the roots with nroots or otherwise you could use numpy's roots function.

In [157]: eqzp = eqz.subs({K1:0.2, K2:0.5})                                                                                                                                                                 

In [158]: eqzp                                                                                                                                                                                              
Out[158]: 
 3                  2                                              
z ⋅(0.6 + 4.5⋅ⅈ) + z ⋅(0.2 + 0.5⋅ⅈ) + z⋅(0.2 - 0.5⋅ⅈ) + 0.6 - 4.5⋅ⅈ

In [159]: Poly(eqzp, z).nroots()                                                                                                                                                                            
Out[159]: [-0.617215947987055 + 0.786793793538333⋅ⅈ, -0.491339121039621 - 0.870968350823388⋅ⅈ, 0.993562347047054 + 0.113286638798883⋅ⅈ]

In [163]: coeffs = [complex(c) for c in Poly(eqzp, z).all_coeffs()]                                                                                                                                            

In [164]: np.roots(coeffs)                                                                                                                                                                                  
Out[164]: 
array([ 0.99356235+0.11328664j, -0.61721595+0.78679379j,
       -0.49133912-0.87096835j])

Either way this gives you 3 possible values for z which is exp(I*theta) so you can get theta (modulo 2*pi) with:

In [167]: r1, r2, r3 = Poly(eqzp, z).nroots()                                                                                                                                                               

In [168]: get_theta = lambda r: acos((r + r.conjugate())/2)                                                                                                                                                 

In [169]: get_theta(r1)                                                                                                                                                                                     
Out[169]: 2.23599562043958

In [170]: get_theta(r2)                                                                                                                                                                                     
Out[170]: 2.08442292239622

In [171]: get_theta(r3)                                                                                                                                                                                     
Out[171]: 0.113530366549989

The transformations we've done mean that +- these values can be solutions to the original equation so we can check by substituting in e.g.:

In [178]: eq.subs({K1:0.2, K2:0.5}).subs(theta, get_theta(r1))                                                                                                                                              
Out[178]: -5.55111512312578e-17

In [179]: eq.subs({K1:0.2, K2:0.5}).subs(theta, get_theta(r2))                                                                                                                                              
Out[179]: -0.124767626702216

In [180]: eq.subs({K1:0.2, K2:0.5}).subs(theta, -get_theta(r2))                                                                                                                                             
Out[180]: 5.55111512312578e-17
user_94509910
Related