【发布时间】:2020-09-10 07:49:53
【问题描述】:
我有下面的工作功能。我有一个函数,我可以从中计算一阶和二阶导数。然后我需要找到一阶导数为零且二阶导数为负的 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)
【问题讨论】:
-
我建议看看
multiprocessing module -
我也想过这个问题,但希望有人能告诉我一个更有效的方法来解决这个问题。如果点数真的很高,多处理只会让我到目前为止。