在 Python 中绘制 Voigt 函数

计算科学 Python scipy 一体化
2021-11-26 21:26:11

我一直在尝试在 Python 中绘制以下函数：

H(a,u)=aπ∫−∞∞exp(−y2)a2+(u−y)2dy

但我不断收到以下错误：

File "//anaconda/lib/python3.5/site-packages/scipy/integrate/quadpack.py", line 383, in _quad
return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

error: Supplied function does not return a valid float.

这是我正在使用的代码：

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

a = 0.1
pi = np.pi
u = np.linspace(-100,100,200)


def integrand(a,u,y):
    top = np.exp(-y**2)
    bottom = a**2 + (u - y)**2
    return top/bottom

def func(u):
    res = np.zeros_like(u)
    for i,val in enumerate(u):
        y,err = integrate.quad(integrand, -np.inf, np.inf, args=(a,u))
        res[i] = y
    return (a/pi)*res


plt.plot(u,func(u))

u 是代表频率范围的一系列值。我一直在努力让集成在 python 中工作。有没有人可以推荐我去让它工作的好方向？

谢谢

2个回答

你看过 的文档scipy.integrate.quad吗？我发现您的代码至少有两个问题：

  • 定义中的第一个参数integrand应该是积分变量y。
  • uin的参数integrand是一个标量，但您在 的args选项中提供了一个向量integrate.quad。

如果您只想绘制函数，则以下代码段简洁明了

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

def Voigt(a, u):
    I = integrate.quad(lambda y: np.exp(-y**2)/(a**2 + (u - y)**2),-np.inf, np.inf)[0]

    return (a/np.pi)*I

a = 0.1
u_range = np.linspace(-100,100,200)

plt.plot(u_range, [Voigt(a, u) for u in u_range])

在此处输入图像描述

您的代码大部分是正确的。只需更改u为val您致电时integrate.quad：

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

a = 0.1
pi = np.pi
u = np.linspace(-100,100,200)


def integrand(a,u,y):
    top = np.exp(-y**2)
    bottom = a**2 + (u - y)**2
    return top/bottom

def func(u):
    res = np.zeros_like(u)
    for i,val in enumerate(u):
        y,err = integrate.quad(integrand, -np.inf, np.inf, args=(a,val))
        res[i] = y
    return (a/pi)*res

plt.plot(u,func(u))

正如@Stelios 所提到的，将u（一个数组）传递给integrand意味着它将返回一个数组，而integrate.quad需要一个标量。在这里，我们只是利用您已经val使用 enumerate.

其它你可能感兴趣的问题
上一篇的定义是什么大号∞L∞- 不合格有限元的规范？ 下一篇使用低等级属性进行最大值/最小值搜索（或排序）