Построение интегрального решения с переменным параметром

Question

Я знаю, что здесь есть похожие вопросы, но ни один из них не подходит к root моей проблемы.

У меня есть интеграл, который содержит параметр, против которого я хочу построить график.

Мой код

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

, если это поможет.

Примечание: mu = x в моем python коде.

Answer 1

Если вы хотите избежать явного l oop, вы можете использовать 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()

Выход: