Python：如何在数字上集成具有两个未知参数的函数

问题描述

现在我分别有两个功能

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）

然后我要针对'u'集成'rho（u）* psi（w（x-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,500)

u = np.linspace(0,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循环的方法。

解决方法

这里有可能使用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,101)
res = empty((x.size,w.size))
for i,xVal in enumerate(x):
    res[i],err = quad_vec(partial(func,x=xVal),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,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,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')