Scipy 的 solve_bvp 和耦合微分方程的性能问题
Performance issue with Scipy's solve_bvp and coupled differential equations
我在 Python 3.8.3 中尝试实现下面的耦合微分方程(也称为单模耦合方程)时遇到问题。至于求解器,我使用的是Scipy的函数scipy.integrate.solve_bvp
,其文档可以阅读here。我想求解复杂域中的方程,对于不同的传播轴值 (z
) 和不同的 beta 值 (beta_analysis
)。
问题在于,与在 Matlab 中使用函数 bvp4c
、bvpinit
和 bvpset
的等效实现相比,它非常慢(不可管理)。评估两个执行的前几次迭代,它们 return 相同的结果,除了在 Scipy 的情况下生成的网格要大得多。网格有时甚至饱和到最大值。
下面显示了要求解的方程以及边界条件函数。
import h5py
import numpy as np
from scipy import integrate
def coupling_equation(z_mesh, a):
ka_z = k # Global
z_a = z # Global
a_p = np.empty_like(a).astype(complex)
for idx, z_i in enumerate(z_mesh):
beta_zf_i = np.interp(z_i, z_a, beta_zf) # Get beta at the desired point of the mesh
ka_z_i = np.interp(z_i, z_a, ka_z) # Get ka at the desired point of the mesh
coupling_matrix = np.empty((2, 2), complex)
coupling_matrix[0] = [-1j * beta_zf_i, ka_z_i]
coupling_matrix[1] = [ka_z_i, 1j * beta_zf_i]
a_p[:, idx] = np.matmul(coupling_matrix, a[:, idx]) # Solve the coupling matrix
return a_p
def boundary_conditions(a_a, a_b):
return np.hstack(((a_a[0]-1), a_b[1]))
此外,鉴于 [=25],我找不到将 k
、z
和 beta_zf
作为函数 coupling_equation
的参数传递的方法=] solve_bpv
函数的参数必须是可调用参数 (x, y)
。我的方法是定义一些全局变量,但如果有更好的解决方案,我也将不胜感激。
我尝试编写的分析函数是:
def analysis(k, z, beta_analysis, max_mesh):
s11_analysis = np.empty_like(beta_analysis, dtype=complex)
s21_analysis = np.empty_like(beta_analysis, dtype=complex)
initial_mesh = np.linspace(z[0], z[-1], 10) # Initial mesh of 10 samples along L
mesh = initial_mesh
# a_init must be complex in order to solve the problem in a complex domain
a_init = np.vstack((np.ones(np.size(initial_mesh)).astype(complex),
np.zeros(np.size(initial_mesh)).astype(complex)))
for idx, beta in enumerate(beta_analysis):
print(f"Iteration {idx}: beta_analysis = {beta}")
global beta_zf
beta_zf = beta * np.ones(len(z)) # Global variable so as to use it in coupling_equation(x, y)
a = integrate.solve_bvp(fun=coupling_equation,
bc=boundary_conditions,
x=mesh,
y=a_init,
max_nodes=max_mesh,
verbose=1)
# mesh = a.x # Mesh for the next iteration
# a_init = a.y # Initial guess for the next iteration, corresponding to the current solution
s11_analysis[idx] = a.y[1][0]
s21_analysis[idx] = a.y[0][-1]
return s11_analysis, s21_analysis
我怀疑问题与传递给不同迭代的初始猜测有关(请参阅 analysis
函数中循环内的注释行)。我尝试将迭代的解决方案设置为以下的初始猜测(这必须减少求解器所需的时间),但它甚至更慢,我不明白。也许我错过了什么,因为这是我第一次尝试求解微分方程。
执行使用的参数如下:
f2 = h5py.File(r'path/to/file', 'r')
k = np.array(f2['k']).squeeze()
z = np.array(f2['z']).squeeze()
f2.close()
analysis_points = 501
max_mesh = 1e6
beta_0 = 3e2;
beta_low = 0; # Lower value of the frequency for the analysis
beta_up = beta_0; # Upper value of the frequency for the analysis
beta_analysis = np.linspace(beta_low, beta_up, analysis_points);
s11_analysis, s21_analysis = analysis(k, z, beta_analysis, max_mesh)
关于如何提高这些功能的性能有什么想法吗?提前谢谢大家,如果问题表述不当,我深表歉意,我接受有关此问题的任何建议。
编辑: 添加了一些有关性能和问题大小的信息。
- 实际上,我找不到确定调用
coupling_equation
次数的关系。这一定是求解器内部操作的问题。我通过打印一行来检查一次迭代中的调用次数,它发生了 133 次(这是最快的一次)。这必须乘以 beta 的迭代次数。对于分析的一个,求解器 return 编辑了这个:
Solved in 11 iterations, number of nodes 529.
Maximum relative residual: 9.99e-04
Maximum boundary residual: 0.00e+00
a
和 z_mesh
的形状是相关的,因为 z_mesh 是一个向量,其长度对应于网格的大小,求解器每次调用 coupling_equation
。鉴于a
包含z_mesh
各点的前进波和回归波的振幅,a
的形状为(2, len(z_mesh))
。
- 在计算时间方面,我只用 Python 在大约 2 小时内完成了 19 次迭代。在这种情况下,初始迭代速度更快,但随着网格的增长,它们开始花费更多时间,直到网格饱和到最大允许值。我认为这是因为那一点的输入耦合系数的值,因为当
beta_analysis
中没有循环被执行时也会发生这种情况(只是 beta 中间值的 solve_bvp
函数)。相反,Matlab 在大约 6 分钟内设法 return 解决了整个问题。如果我将最后一次迭代的结果作为 initial_guess
传递(analysis
函数中的注释行,网格溢出得更快并且不可能获得超过几次迭代。
基于半随机输入,我们可以看到有时会达到 max_mesh
。这意味着可以使用相当大的 z_mesh
和 a
数组调用 coupling_equation
。问题是 coupling_equation
包含一个 slow pure-Python loop 在数组的每一列上迭代。您可以使用 Numpy 向量化 大大加快计算速度。这是一个实现:
def coupling_equation_fast(z_mesh, a):
ka_z = k # Global
z_a = z # Global
a_p = np.empty(a.shape, dtype=np.complex128)
beta_zf_i = np.interp(z_mesh, z_a, beta_zf) # Get beta at the desired point of the mesh
ka_z_i = np.interp(z_mesh, z_a, ka_z) # Get ka at the desired point of the mesh
# Fast manual matrix multiplication
a_p[0] = (-1j * beta_zf_i) * a[0] + ka_z_i * a[1]
a_p[1] = ka_z_i * a[0] + (1j * beta_zf_i) * a[1]
return a_p
与原始实现相比,此代码使用半随机输入提供了类似的输出,但在我的机器上快了大约 20 倍。
此外,我不知道 max_mesh
是否恰好对您的输入也很大,即使这是 normal/intended。减少 max_mesh
的值以进一步减少执行时间可能是有意义的。
我在 Python 3.8.3 中尝试实现下面的耦合微分方程(也称为单模耦合方程)时遇到问题。至于求解器,我使用的是Scipy的函数scipy.integrate.solve_bvp
,其文档可以阅读here。我想求解复杂域中的方程,对于不同的传播轴值 (z
) 和不同的 beta 值 (beta_analysis
)。
问题在于,与在 Matlab 中使用函数 bvp4c
、bvpinit
和 bvpset
的等效实现相比,它非常慢(不可管理)。评估两个执行的前几次迭代,它们 return 相同的结果,除了在 Scipy 的情况下生成的网格要大得多。网格有时甚至饱和到最大值。
下面显示了要求解的方程以及边界条件函数。
import h5py
import numpy as np
from scipy import integrate
def coupling_equation(z_mesh, a):
ka_z = k # Global
z_a = z # Global
a_p = np.empty_like(a).astype(complex)
for idx, z_i in enumerate(z_mesh):
beta_zf_i = np.interp(z_i, z_a, beta_zf) # Get beta at the desired point of the mesh
ka_z_i = np.interp(z_i, z_a, ka_z) # Get ka at the desired point of the mesh
coupling_matrix = np.empty((2, 2), complex)
coupling_matrix[0] = [-1j * beta_zf_i, ka_z_i]
coupling_matrix[1] = [ka_z_i, 1j * beta_zf_i]
a_p[:, idx] = np.matmul(coupling_matrix, a[:, idx]) # Solve the coupling matrix
return a_p
def boundary_conditions(a_a, a_b):
return np.hstack(((a_a[0]-1), a_b[1]))
此外,鉴于 [=25],我找不到将 k
、z
和 beta_zf
作为函数 coupling_equation
的参数传递的方法=] solve_bpv
函数的参数必须是可调用参数 (x, y)
。我的方法是定义一些全局变量,但如果有更好的解决方案,我也将不胜感激。
我尝试编写的分析函数是:
def analysis(k, z, beta_analysis, max_mesh):
s11_analysis = np.empty_like(beta_analysis, dtype=complex)
s21_analysis = np.empty_like(beta_analysis, dtype=complex)
initial_mesh = np.linspace(z[0], z[-1], 10) # Initial mesh of 10 samples along L
mesh = initial_mesh
# a_init must be complex in order to solve the problem in a complex domain
a_init = np.vstack((np.ones(np.size(initial_mesh)).astype(complex),
np.zeros(np.size(initial_mesh)).astype(complex)))
for idx, beta in enumerate(beta_analysis):
print(f"Iteration {idx}: beta_analysis = {beta}")
global beta_zf
beta_zf = beta * np.ones(len(z)) # Global variable so as to use it in coupling_equation(x, y)
a = integrate.solve_bvp(fun=coupling_equation,
bc=boundary_conditions,
x=mesh,
y=a_init,
max_nodes=max_mesh,
verbose=1)
# mesh = a.x # Mesh for the next iteration
# a_init = a.y # Initial guess for the next iteration, corresponding to the current solution
s11_analysis[idx] = a.y[1][0]
s21_analysis[idx] = a.y[0][-1]
return s11_analysis, s21_analysis
我怀疑问题与传递给不同迭代的初始猜测有关(请参阅 analysis
函数中循环内的注释行)。我尝试将迭代的解决方案设置为以下的初始猜测(这必须减少求解器所需的时间),但它甚至更慢,我不明白。也许我错过了什么,因为这是我第一次尝试求解微分方程。
执行使用的参数如下:
f2 = h5py.File(r'path/to/file', 'r')
k = np.array(f2['k']).squeeze()
z = np.array(f2['z']).squeeze()
f2.close()
analysis_points = 501
max_mesh = 1e6
beta_0 = 3e2;
beta_low = 0; # Lower value of the frequency for the analysis
beta_up = beta_0; # Upper value of the frequency for the analysis
beta_analysis = np.linspace(beta_low, beta_up, analysis_points);
s11_analysis, s21_analysis = analysis(k, z, beta_analysis, max_mesh)
关于如何提高这些功能的性能有什么想法吗?提前谢谢大家,如果问题表述不当,我深表歉意,我接受有关此问题的任何建议。
编辑: 添加了一些有关性能和问题大小的信息。
- 实际上,我找不到确定调用
coupling_equation
次数的关系。这一定是求解器内部操作的问题。我通过打印一行来检查一次迭代中的调用次数,它发生了 133 次(这是最快的一次)。这必须乘以 beta 的迭代次数。对于分析的一个,求解器 return 编辑了这个:
Solved in 11 iterations, number of nodes 529. Maximum relative residual: 9.99e-04 Maximum boundary residual: 0.00e+00
a
和z_mesh
的形状是相关的,因为 z_mesh 是一个向量,其长度对应于网格的大小,求解器每次调用coupling_equation
。鉴于a
包含z_mesh
各点的前进波和回归波的振幅,a
的形状为(2, len(z_mesh))
。- 在计算时间方面,我只用 Python 在大约 2 小时内完成了 19 次迭代。在这种情况下,初始迭代速度更快,但随着网格的增长,它们开始花费更多时间,直到网格饱和到最大允许值。我认为这是因为那一点的输入耦合系数的值,因为当
beta_analysis
中没有循环被执行时也会发生这种情况(只是 beta 中间值的solve_bvp
函数)。相反,Matlab 在大约 6 分钟内设法 return 解决了整个问题。如果我将最后一次迭代的结果作为initial_guess
传递(analysis
函数中的注释行,网格溢出得更快并且不可能获得超过几次迭代。
基于半随机输入,我们可以看到有时会达到 max_mesh
。这意味着可以使用相当大的 z_mesh
和 a
数组调用 coupling_equation
。问题是 coupling_equation
包含一个 slow pure-Python loop 在数组的每一列上迭代。您可以使用 Numpy 向量化 大大加快计算速度。这是一个实现:
def coupling_equation_fast(z_mesh, a):
ka_z = k # Global
z_a = z # Global
a_p = np.empty(a.shape, dtype=np.complex128)
beta_zf_i = np.interp(z_mesh, z_a, beta_zf) # Get beta at the desired point of the mesh
ka_z_i = np.interp(z_mesh, z_a, ka_z) # Get ka at the desired point of the mesh
# Fast manual matrix multiplication
a_p[0] = (-1j * beta_zf_i) * a[0] + ka_z_i * a[1]
a_p[1] = ka_z_i * a[0] + (1j * beta_zf_i) * a[1]
return a_p
与原始实现相比,此代码使用半随机输入提供了类似的输出,但在我的机器上快了大约 20 倍。
此外,我不知道 max_mesh
是否恰好对您的输入也很大,即使这是 normal/intended。减少 max_mesh
的值以进一步减少执行时间可能是有意义的。