我的优化任务涉及以下积分的计算，并找到xl和xu的最佳值：
迭代花费的时间太长，因此我决定通过为所有可能的值xl和xu计算积分来加快它们的速度，然后在优化过程中对计算出的值进行插值。
我写了以下函数：
def k_integrand(x, xl, xu):
return((x**2)*mpmath.exp(x))/((xu - xl)*(mpmath.exp(x)-1)**2)
@np.vectorize
def K(xl, xu):
y, err = integrate.quad(k_integrand, xl, xu, args = (xl, xu))
return y
和两个相同的数组
grid_xl和grid_xu，其值以动态递增。当我运行代码时，我得到以下信息：
K(grid_xl, grid_xu)
Traceback (most recent call last):
File "<ipython-input-2-5b9df02f12b7>", line 1, in <module>
K(grid_xl, grid_xu)
File "C:/Users/909404/OneDrive/Работа/ZnS-FeS/Теплоемкость/Python/CD357/4 - Optimization CD357 interpolation.py", line 75, in K
y, err = integrate.quad(k_integrand, xl, xu, args = (xl, xu))
File "C:\Users\909404\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py", line 323, in quad
points)
File "C:\Users\909404\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py", line 372, in _quad
if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
我想这是因为
xl应该总是小于xu。有什么方法可以比较
xl和xu的值并在xl>=xu的情况下返回NaN？最后，我想要这样的东西：
并具有使用插值的能力。
也许我选择了错误的方式？我将不胜感激。
最佳答案
除非省略np.vectorize装饰器，否则无法重现您的错误。设置一致的xl / xu值确实会给我一个ZeroDivisionError。
无论如何，没有什么可以阻止您检查更高级别函数中的xu与xl的值。这样，您可以完全跳过对无意义的数据点的集成，并尽早返回np.nan：
import numpy as np
import mpmath
import scipy.integrate as integrate
def k_integrand(x, xl, xu):
return ((x**2)*mpmath.exp(x))/((xu - xl)*(mpmath.exp(x)-1)**2)
@np.vectorize
def K(xl, xu):
if xu <= xl:
# don't even try to integrate
return np.nan
y, err = integrate.quad(k_integrand, xl, xu, args = (xl, xu))
return y
grid_xl = np.linspace(0.1,1,10) # shape (10,) ~ (1,10)
grid_xu = np.linspace(0.5,4,8)[:,None] # shape (8,1)
有了这些定义，我得到了（以下
np.set_printoptions(linewidth=200)为便于比较：In [35]: K(grid_xl, grid_xu)
Out[35]:
array([[0.99145351, 0.98925197, 0.98650808, 0.98322919, nan, nan, nan, nan, nan, nan],
[0.97006703, 0.96656815, 0.96254363, 0.95800307, 0.95295785, 0.94742104, 0.94140733, 0.93493293, 0.9280154 , nan],
[0.93730403, 0.93263063, 0.92745487, 0.92178832, 0.91564423, 0.90903747, 0.90198439, 0.89450271, 0.88661141, 0.87833062],
[0.89565597, 0.88996696, 0.88380385, 0.87717991, 0.87010995, 0.8626103 , 0.85469862, 0.84639383, 0.83771595, 0.82868601],
[0.84794429, 0.8414176 , 0.83444842, 0.82705134, 0.81924245, 0.81103915, 0.8024601 , 0.79352503, 0.7842547 , 0.77467065],
[0.79692339, 0.78974 , 0.78214742, 0.77416128, 0.76579857, 0.75707746, 0.74801726, 0.73863822, 0.72896144, 0.71900874],
[0.7449893 , 0.73732055, 0.7292762 , 0.72087263, 0.71212741, 0.70305921, 0.69368768, 0.68403329, 0.67411725, 0.66396132],
[0.69402415, 0.68602325, 0.67767956, 0.66900991, 0.66003222, 0.65076537, 0.6412291 , 0.63144388, 0.62143077, 0.61121128]])
您可以看到这些值与链接的图像完全一致。
现在，我有一个坏消息和一个好消息。坏消息是，尽管
np.vectorize提供了围绕数组输入调用标量集成函数的语法糖，但与本机for循环相比，它实际上并没有提高速度。好消息是，您可以将对mpmath.exp的调用替换为对np.exp的调用，并且最终将获得相同的结果：def k_integrand_np(x, xl, xu):
return ((x**2)*np.exp(x))/((xu - xl)*(np.exp(x)-1)**2)
@np.vectorize
def K_np(xl, xu):
if xu <= xl:
# don't even try to integrate
return np.nan
y, err = integrate.quad(k_integrand_np, xl, xu, args = (xl, xu))
return y
有了这些定义
In [14]: res_mpmath = K(grid_xl, grid_xu)
...: res_np = K_np(grid_xl, grid_xu)
...: inds = ~np.isnan(res_mpmath)
...:
In [15]: np.array_equal(res_mpmath[inds], res_np[inds])
Out[15]: True
In [16]: %timeit K(grid_xl, grid_xu)
107 ms ± 521 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [17]: %timeit K_np(grid_xl, grid_xu)
7.26 ms ± 157 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
因此，这两种方法给出的结果相同（完全是！），但是numpy版本的速度快将近15倍。
关于python - 计算2D插值的积分时出错。比较numpy数组，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/49557809/