Python绘制Airy函数

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

y ′ ′ ( x ) = x y ( 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函数的积分形式可写为

A i ( x ) = 1 π ∫ cos ⁡ ( 1 3 t 3 + x t ) d t B i ( x ) = 1 π ∫ e − t 3 3 + x t + sin ⁡ ( t 3 3 + x t ) d t begin{aligned} A_i(x)&=frac{1}{pi}intcos(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(31​t3+xt)dt=π1​∫e−3t3​+xt+sin(3t3​+xt)dt​

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

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

其中，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、Bi及其导数的零点所在的位置。以ai_zeros(n)为例，其返回值共有四组

• a： A i ( x ) A_i(x) Ai​(x)的前n个零点的位置
• ap： A i ′ ( x ) A_i'(x) Ai′​(x)的前n个零点的位置
• ai： A i ( x ) A_i(x) Ai​(x)的前n个零点对应的横坐标在 A i ′ ( x ) A_i'(x) Ai′​(x)上的值
• aip： A i ′ ( x ) A_i'(x) Ai′​(x)的前n个零点对应的横坐标在 A i ( 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(x)为Airy函数的积分，返回值分别是

A p t = ∫ 0 x A i ( t ) d t B p t = ∫ 0 x B i ( t ) d t A n t = ∫ 0 x A i ( − t ) d t B n t = ∫ 0 x A i ( − t ) d t begin{aligned} A_{pt}=int_0^xA_i(t)text dtquad&B_{pt}=int_0^xB_i(t)text dt\ A_{nt}=int_0^xA_i(-t)text dtquad&B_{nt}=int_0^xA_i(-t)text dt end{aligned} Apt​=∫0x​Ai​(t)dtAnt​=∫0x​Ai​(−t)dt​Bpt​=∫0x​Bi​(t)dtBnt​=∫0x​Ai​(−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()