绘制变参数积分解

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轴上画成一个变量。我该怎么做

在mathematica中，我能做到这一点，情节是这样的

我的mathematica代码是

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}]

如果有帮助的话

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

Answer 1

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()

输出：

Answer 2

如果希望避免显式循环，可以使用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()