如何解决用fsolve解决数组或函数列表的最快方法
我具有以下工作功能。我有一个函数可以从中计算一阶和二阶导数。然后，我需要找到theta的值，其中一阶导数为零，第二个导数为负。我必须针对大量点进行计算。点的数量等于K1和K2的长度。我使用sympy计算一阶和二阶导数。我目前遍历所有导数，并为每个导数求解方程。有没有一种更快的方法，一旦K1和K2的长度增加> 1000，这对我的应用来说就花了很多时间。
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
输入示例：
K1 = np.random.rand(10000,)
K2 = np.random.rand(10000,)
get_cpd(K1,K2)
解决方法
最好的办法是尝试根据符号参数对符号进行尽可能多的符号处理。有可能获得例如first_derivative，但您需要对其进行一些转换。在这里，我将sin / cos重写为exp，然后使用替换exp(I*theta/2) = sqrt(z)来获得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₂)
现在sympy可以解决这个问题（roots(eqz,z)），但是三次方程的通用公式非常复杂，因此可能不是最佳方法。给定K1和K2的特殊浮点值，尽管sympy可以很容易地以nroots求根，否则您可以使用numpy的roots函数。
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])
这两种方式都可以为z提供3个可能的值，即exp(I*theta)，因此您可以通过以下方式获得theta（模2*pi）：
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
我们完成的转换意味着+-这些值可以作为原始方程的解决方案，因此我们可以通过以下方式代替：
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,get_theta(r2))
Out[179]: -0.124767626702216
In [180]: eq.subs({K1:0.2,-get_theta(r2))
Out[180]: 5.55111512312578e-17
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。