Python：如何对具有两个未知参数的函数进行数值积分

Question

现在我有两个功能分别是

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

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

然后我想整合 'rho(u) * psi( w(x-u) )' 关于 'u'。使得积分结果可以是关于 'w' 和 'x'.

的一个函数

这是我尝试求解此积分时的 Python 代码片段。

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 中的 multi-dimensional 数组。下面更新的方法可以将 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')

Python：如何对具有两个未知参数的函数进行数值积分

Python: how to integrate functions with two unknown parameters numerically

python

numpy

scipy

numerical-integration