【Python】scipy之艾里函数

文章目录

Airy函数
积分形式
Airy相关函数
- airye
- ai_zeros，bi_zeros
- itairy

Airy函数

Ai和Bi是微分方程的两个独立的解

y′′(x)=xy(x)y''(x)=xy(x) y′′(x)=xy(x)

函数airy的返回值有4个，分别是Ai, Bi及其导数。

import scipy.special as sc
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(-15,2,0.1)
Ai, Aip, Bi, Bip = sc.airy(x)
plt.plot(x,Ai)
plt.plot(x,Bi)
plt.show()

如图所示

积分形式

Airy函数的积分形式可写为

Ai(x)=1π∫cos⁡(13t3+xt)dtBi(x)=1π∫e−t33+xt+sin⁡(t33+xt)dt\begin{aligned} A_i(x)&=\frac{1}{\pi}\int\cos(\frac{1}{3}t^3+xt)\text dt\\ B_i(x)&=\frac{1}{\pi}\int e^{-\frac{t^3}{3}+xt}+\sin(\frac{t^3}{3}+xt) \text dt \end{aligned} Ai(x)Bi(x)=π1∫cos(31t3+xt)dt=π1∫e−3t3+xt+sin(3t3+xt)dt

Airy相关函数

在scipy的特殊函数包中，除了airy之外，还封装了四个与Airy函数相关的函数

airye：指数缩放的Airy函数及其导数
ai_zeros：Ai及其导数的零点
bi_zeros：Bi及其导数的零点
itairy：Airy函数的积分

airye

其中，airye与airy的缩放关系如下

k = 2.0/3.0
eAi  = Ai  * np.exp(k*z*sqrt(z))
eAip = Aip * np.exp(k*z*sqrt(z))
eBi  = Bi  * np.exp(-abs(k*(z*sqrt(z)).real))
eBip = Bip * np.exp(-abs(k*(z*sqrt(z)).real))

绘制图像如下

代码为

z = np.linspace(0, 50, 500)
data = sc.airye(z)
f, axes = plt.subplots(2, 1)
names = [["eAi", "eAip"], ["eBi", "eBip"]]
for i in range(2):axes[i].plot(z, data[2*i], "-r", z, data[2*i+1], "-b")axes[i].legend(names[i])axes[i].grid(True)plt.show()

ai_zeros，bi_zeros

ai_zeros和bi_zeros输入为零点个数，输出为Ai、Bi及其导数的零点所在的位置。以ai_zeros(n)为例，其返回值共有四组

a：Ai(x)A_i(x)Ai(x)的前n个零点的位置
ap：Ai′(x)A_i'(x)Ai′(x)的前n个零点的位置
ai：Ai(x)A_i(x)Ai(x)的前n个零点对应的横坐标在Ai′(x)A_i'(x)Ai′(x)上的值
aip：Ai′(x)A_i'(x)Ai′(x)的前n个零点对应的横坐标在Ai(x)A_i(x)Ai(x)上的值

这两个函数与airy的关系如下

绘图代码为

x = np.arange(-15,2,0.1)
data = sc.airy(x)
Ai, Aip, Bi, Bip = sc.airy(x)
names = [["eAi", "eAip"], ["eBi", "eBip"]]
f, axes = plt.subplots(2, 1)
# 绘制Ai和Bi及其导数的图像
for i in range(2):axes[i].plot(x, data[2*i], "-r", x, data[2*i+1], "-b")axes[i].legend(names[i])axes[i].grid(True)# 求出Airy函数的相应零点
aZeros = sc.ai_zeros(5)
bZeros = sc.bi_zeros(5)axes[0].stem(aZeros[0],aZeros[3])
axes[0].stem(aZeros[1],aZeros[2])
axes[1].stem(bZeros[0],bZeros[3])
axes[1].stem(bZeros[1],bZeros[2])plt.show()

itairy

最后，itairy(x)为Airy函数的积分，返回值分别是

Apt=∫0xAi(t)dtBpt=∫0xBi(t)dtAnt=∫0xAi(−t)dtBnt=∫0xAi(−t)dt\begin{aligned} A_{pt}=\int_0^xA_i(t)\text dt\quad&B_{pt}=\int_0^xB_i(t)\text dt\\ A_{nt}=\int_0^xA_i(-t)\text dt\quad&B_{nt}=\int_0^xA_i(-t)\text dt \end{aligned} Apt=∫0xAi(t)dtAnt=∫0xAi(−t)dtBpt=∫0xBi(t)dtBnt=∫0xAi(−t)dt

其前两个返回值的图像为

绘图代码为

x = np.arange(-15,2,0.1)
data = sc.itairy(x)
plt.plot(x, data[0], "-r", x, data[1], "-b")
plt.legend(["Apt", "Bpt"])
plt.grid(True)plt.show()