绘制具有变化参数的积分解

Question

我知道这里也有类似的问题，但没有一个能触及我问题的根源。

我有一个积分，其中包含一个参数，我想根据该参数生成绘图。

我的代码是

import numpy as np

import matplotlib.pyplot as plt

from scipy import integrate


def intyfun(x, a):
    return np.exp(-a/4)*np.exp(-x**2*a)*(2*np.pi*x)*(np.sinh(np.pi)/np.cosh(np.pi*x)**2)

现在我被困住了。我想对 x 在 0 到无穷大的范围内积分该函数，并将其值绘制为 x 轴上的变化作为参数。我该怎么做？

在数学中我可以做到这一点，情节看起来像这样

我的数学代码是

integral[a_?NumericQ] := 
 NIntegrate[
  Exp[-a/4]*Exp[-mu^2*a]*(2*Pi*mu*Sinh[mu*Pi])/(Cosh[mu*Pi]^2), {mu, 
   0, Infinity}, 
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
    "MaxErrorIncreases" -> 10000, "SingularityHandler" -> "IMT"}, 
  MaxRecursion -> 100, PrecisionGoal -> 4]

Plot[integral[a], {a, 0.01, 10}, ImageSize -> Large, 
 FrameStyle -> Black,
 BaseStyle -> {FontFamily -> "Latin Modern Roman"}, PlotLabel -> "", 
 PlotStyle -> Black, FrameStyle -> Black,
 BaseStyle -> {FontFamily -> "Latin Modern Roman"}, PlotRange -> All, 
 AxesLabel -> {a, IntegralValue}]

如果有帮助的话。

N.B mu=x 在我的 python 代码中。

Answer 1

如果您想避免显式循环，可以使用 quadpy（我的一个项目）在一个矢量化积分步骤中计算所有值。这快得多：

import quadpy
import numpy as np
import matplotlib.pyplot as plt


a = np.linspace(0.0, 10.0, 300)


def intyfun(x):
    return (
        np.exp(-a / 4)[:, None]
        * np.exp(np.multiply.outer(a, -(x ** 2)))
        * (2 * np.pi * x)
        * (np.sinh(np.pi) / np.cosh(np.pi * x) ** 2)
    )


val, _ = quadpy.quad(intyfun, 0, np.inf)

plt.plot(a, val)
plt.grid()
plt.gca().set_aspect("equal")
plt.show()

Answer 2

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


def function(x, a):
    return np.exp(-a/4)*np.exp(-x**2*a)*(2*np.pi*x)*(np.sinh(np.pi)/np.cosh(np.pi*x)**2)


def integrate_function(x_min, x_max, a):
    return integrate.quad(function, x_min, x_max, args=(a,))[0]


# define integration range
x_min = 0
x_max = 1
# define range of a variable
a_points = np.linspace(0, 10, 100)

results = []
for a in a_points:
    value = integrate_function(x_min, x_max, a)
    results.append(value)


plt.plot(a_points, results)
plt.ylabel("Integral value")
plt.xlabel("a variable")
plt.show()

输出：