Python：如何以数值方式集成具有两个未知参数的函数

Question

现在我有两个函数分别是

rho(u) = np.exp((-2.0 / 0.2) * (u**0.2-1.0))

psi( w(xu) ) = (1/(4.0 * math.sqrt(np.pi))) * np.exp(- ((w * (xu))**2) / 4.0) * ( 2.0 - (w * (xu))**2)

然后我想将 'rho(u) * psi( w(x-u) )' 与 'u' 积分。这样积分结果可以是关于'w'和'x'的一个函数。

这是我尝试求解这个积分时的 Python 代码 sn-p。

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import integrate

x = np.linspace(0,10,1000)
w = np.linspace(0,10,500)

u = np.linspace(0,10,1000)

rho = np.exp((-2.0/0.2)*(u**0.2-1.0))

value = np.zeros((500,1000),dtype="float32")

# Integrate the products of rho with 
# (1/(4.0*math.sqrt(np.pi)))*np.exp(- ((w[i]*(x[j]-u))**2) / 4.0)*(2.0 - (w[i]*(x[j]-u))**2)
for i in range(len(w)):
    for j in range(len(x)):
        value[i,j] =value[i,j]+ integrate.simps(rho*(1/(4.0*math.sqrt(np.pi)))*np.exp(- ((w[i]*(x[j]-u))**2) / 4.0)*(2.0 - (w[i]*(x[j]-u))**2),u)

plt.imshow(value,origin='lower')
plt.colorbar()

如上所示，当我进行集成时，我使用了嵌套 for 循环。我们都知道这样的方式是低效的。

所以想问问有没有不使用for循环的方法。

Answer 1

这是使用scipy.integrate.quad_vec 的可能性。它在我的机器上执行 6 秒，我认为这是可以接受的。不过，我确实只对x 和w 使用了0.1 的步长，但这样的分辨率似乎是单核上的一个很好的折衷方案。

from functools import partial
import matplotlib.pyplot as plt
from numpy import empty, exp, linspace, pi, sqrt
from scipy.integrate import quad_vec
from time import perf_counter


def func(u, x):
    rho = exp(-10 * (u ** 0.2 - 1))
    var = w * (x - u)
    psi = exp(-var ** 2 / 4) * (2 - var ** 2) / 4 / sqrt(pi)
    return rho * psi


begin = perf_counter()
x = linspace(0, 10, 101)
w = linspace(0, 10, 101)
res = empty((x.size, w.size))
for i, xVal in enumerate(x):
    res[i], err = quad_vec(partial(func, x=xVal), 0, 10)
print(f'{perf_counter() - begin} s')
plt.contourf(w, x, res)
plt.colorbar()
plt.xlabel('w')
plt.ylabel('x')
plt.show()

更新

我没有意识到，但也可以使用quad_vec 中的多维数组。下面的更新方法可以将x 和w 的分辨率提高2 倍，并保持大约7 秒的执行时间。此外，没有更多可见的for 循环。

import matplotlib.pyplot as plt
from numpy import exp, mgrid, pi, sqrt
from scipy.integrate import quad_vec
from time import perf_counter


def func(u):
    rho = exp(-10 * (u ** 0.2 - 1))
    var = w * (x - u)
    psi = exp(-var ** 2 / 4) * (2 - var ** 2) / 4 / sqrt(pi)
    return rho * psi


begin = perf_counter()
x, w = mgrid[0:10:201j, 0:10:201j]
res, err = quad_vec(func, 0, 10)
print(f'{perf_counter() - begin} s')
plt.contourf(w, x, res)
plt.colorbar()
plt.xlabel('w')
plt.ylabel('x')
plt.show()

处理评论

只需在plt.show() 之前添加以下行，即可使两个轴以对数方式缩放。

plt.gca().set_xlim(0.05, 10)
plt.gca().set_ylim(0.05, 10)
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')