1 特殊函数 (scipy.special)
本专栏主要按照SciPy官网的Tutorial介绍SciPy的各种子库及其应用。
说起科学计算，肯定少不了各种特殊函数。作为一个物理系学生，我对各种特殊函数可是如数家珍。本篇文章主要介绍一下物理中常用的特殊函数，顺带着复习一下数学物理方法吧，毕竟特殊函数基本上都是从解微分方程解出来的。
本文中使用的特殊函数都需提前导入SciPy的特殊函数库，有些地方还需使用NumPy库与Matplotlib。
import numpy as np
import matplotlib.pyplot as plt
from scipy import special1 艾里函数
求解艾里方程时：
y''=xy\\
可得两个解：
Ai(x)=\frac1\pi\int_{0}^{+\infty}\cos\left(\frac{t^3}{3}+xt\right)\text{d}t\\ Bi(x)=\frac1\pi\int_{0}^{+\infty}\left[e^{-\frac{t^3}{3}+xt}+\sin(\frac{t^3}{3}+xt)\right]\text{d}t\\
Ai(x)称为（第一）艾里函数，Bi(x)称为第二艾里函数。
Scipy中的艾里函数：
special.airy()输入的值可以为一个数：
print(special.airy(1))返回的是一个数组：
(0.13529241631288147, -0.15914744129679328, 1.2074235949528715, 0.9324359333927756)其中分别是Ai(1)，Aip(1)，Bi(1)，Bip(1)的值，Aip、Bip为Ai、Bi的导数。
输入还可为一个数组：
x=np.linspace(0,1,10)
print(special.airy(x))返回的是：
(
#Ai(x)
array([0.35502805, 0.32634823, 0.2981096 , 0.27068265, 0.24436385,
0.21938061, 0.19589735, 0.17402233, 0.15381484, 0.13529242]),
#Aip(x)
array([-0.2588194 , -0.25674602, -0.25099405, -0.24224712, -0.2311511 ,
-0.21830163, -0.20423569, -0.18942675, -0.17428288, -0.15914744]),
#Bi(x)
array([0.61492663, 0.66488273, 0.71576249, 0.76861822, 0.82465049,
0.88522296, 0.9518849 , 1.02640275, 1.11080234, 1.20742359]),
Bip(x)
array([0.44828836, 0.45228954, 0.46512342, 0.48807833, 0.52254465,
0.57007899, 0.63247928, 0.71187269, 0.81081937, 0.93243593])
)其中每个数组的值与Ai，Aip，Bi，Bip相对应。
这样就可据此作图：
x=np.linspace(-10,5,1000)
plt.plot(x,special.airy(x)[0],label='Ai(x)')
plt.plot(x,special.airy(x)[2],label='Bi(x)')
plt.legend()
plt.grid()
plt.ylim(-0.5,1)
plt.show()
除此之外，SciPy还提供了艾里函数的积分：
special.itairy(x) # x can be any real number.该函数返回一个数组：
(Apt, Bpt, Ant, Bnt)其分别代表：
\text{Apt}=\int_{0}^{x}Ai(t)\text{d}t\\ \text{Bpt}=\int_{0}^{x}Bi(t)\text{d}t\\ \text{Ant}=\int_{0}^{x}Ai(-t)\text{d}t\\ \text{Bnt}=\int_{0}^{x}Bi(-t)\text{d}t\\
2 椭圆函数
1.雅可比椭圆函数
这里可能要先介绍一下雅可比椭圆函数。设：
u=\int_{0}^{\phi}\frac{\text{d}\theta}{\sqrt{1-m\sin^2\theta}}\\
椭圆正弦函数（sn u）、椭圆余弦函数（cn u）、椭圆德尔塔函数（dn u）定义为：
\text{sn}\ u=\sin\phi\\ \text{cn}\ u=\cos\phi\\ \text{dn}\ u=\sqrt{1-m\sin^2\phi}
SciPy中的雅可比椭圆函数是：
special.ellipj(u,m) # m is a parameter between 0 and 1.这里m的位置如果输入数组，默认取数组的第一个数。
返回值是：
(sn(u|m),cn(u|m), dn(u|m), phi)其中phi即雅可比椭圆函数定义中的 \phi 。
2.第一类完全椭圆积分
第一类完全椭圆积分的表达式为：
K(k)=\int_{0}^{\frac\pi2}\frac{\text{d}\theta}{\sqrt{1-k\sin^2\theta}}\\
SciPy中提供的函数为：
special.ellipk(k) # k is the independent variable.其中k可为数也可为数组，返回与输入相同类型的数据。
SciPy中还提供了这样的函数：
special.ellipkm1(p) # p is the independent variable.其中p可为数也可为数组，返回与输入相同类型的数据。
这个函数专门用于计算参数k在1附近时的情况。具体来说，它计算的是：
K(p)=\int_{0}^{\frac{\pi}{2}}\frac{\text{d}\theta}{\sqrt{1-(1-p)\sin^2\theta}}\\
3.第一类不完全椭圆积分
第一类不完全椭圆积分的表达式为：
F(\phi|m)=\int_0^\phi\frac{\text{d}\theta}{\sqrt{1-m\sin^2\theta}}\\
SciPy提供的函数为：
special.ellipkinc(phi,m)4.第二类完全椭圆积分
第二类完全椭圆积分的表达式为：
E(k)=\int_0^\frac\pi2\sqrt{1-k\sin^2\theta}\text{d}\theta\\
SciPy提供的函数为：
special.ellipe(k)5.第二类不完全椭圆积分
第二类不完全椭圆积分的表达式为：
E(\phi,k)=\int_0^\phi\sqrt{1-k\sin^2\theta}\text{d}\theta\\
SciPy提供的函数为：
special.ellipeinc(phi,k)3 贝塞尔函数
1.第一、二类贝塞尔函数
贝塞尔函数是学光学时就认识的老朋友了，大名鼎鼎的艾里斑的推导结果中就有贝塞尔函数。
点光源通过理想透镜成像，其衍射光强有相应的数学模型：
I(\theta)=I_0\left(\frac{2J_1(ka\sin\theta)}{ka\sin\theta}\right)^2=I_0\left(\frac{2J_1(x)}{x}\right)^2\\
这里的 J_1(x) 就是第一类贝塞尔函数。
Scipy提供的贝塞尔函数为：
special.jv(v,z)它计算的表达式是：
J_v(z)=\sum_{m=0}^{+\infty}\frac{(-1)^m}{m!\Gamma(m+v+1)}\left(\frac{x}{2}\right)^{2m+v}\\
这里的z可以为复数。
对于艾里斑的情况，我们可以做出 J_1(z) 的图像：
可以发现和实验得到的艾里斑符合的很好，光强基本都集中在中心处，周围迅速变暗。
事实上，贝塞尔函数是下列方程（贝塞尔方程）的解：
x^2y\frac{\text{d}^2y}{\text{d}x^2}+x\frac{\text{d}y}{\text{d}x}+(x^2-v^2)y=0\\
该方程的通解为：
y(x)=c_1J_v(x)+c_2Y_v(x)\\
其中 J_v(x) 称为第一类贝塞尔函数， Y_v(x) 称为第二类贝塞尔函数（或称诺依曼函数），上面计算的是第一类贝塞尔函数。
第二类贝塞尔函数可表示为：
Y_v(z)=\frac{J_v(z)\cos(\pi v)-J_{-v}(z)}{\sin(\pi v)}\\
SciPy提供的函数为：
special.yv(v,z)2.汉克尔函数
汉克尔函数又称第三类贝塞尔函数，其表达式为：
H_v^{(1)}(z)=J_v(z)+iY_v(z)\\ H_v^{(2)}(z)=J_v(z)-iY_v(z)\\
SciPy提供的函数分别为：
special.hankel1(v,z)
special.hankel2(v,z)3.球贝塞尔函数
使用分离变量法求解球坐标下的三维亥姆霍兹方程时，可得到径向方程为：
x^2\frac{\text{d}^2y}{\text{d}x^2}+2x\frac{\text{d}y}{\text{d}x}+[x^2-n(n+1)]y=0\\
其两线性无关的解为：
j_n(x)=\sqrt{\frac{\pi}{2x}}J_{n+\frac12}(x)\\ y_n(x)=\sqrt{\frac{\pi}{2x}}Y_{n+\frac12}(x)\\
SciPy提供的函数为：
special.spherical_jn(n,x)
special.spherical_yn(n,x)4 伽马函数
对于正整数n，伽马函数的表达式为：
\Gamma(n)=(n-1)!\\
将阶乘推广到复数域上可得伽马函数的表达式：
\Gamma(z)=\int_0^{+\infty}t^{z-1}e^{-t}\text{d}t\\
SciPy提供的函数为：
special.gamma(z)当然对于正整数算阶乘的话，还是更推荐直接使用阶乘函数：
special.factorial(n)Scipy还提供了以下函数：
f_1(x)=\ln(|\Gamma(x)|)\\ f_2(x)=\left\{\begin{array} +1,\quad \Gamma(x)>0\\ -1,\quad \Gamma(x)<0\\ \end{array}\right.\\ P(a,x)=\frac{1}{\Gamma(a)}\int_0^{x}t^{\alpha-1}e^{-t}\text{d}t\\ B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\\ f_3(a,b)=\ln(|B(a,b)|)\\
special.gammaln(x) #f_1(x)
special.gammasgn(x) #f_2(x)
special.gammainc(a,x) #P(a,x)
special.beta(a,b) #B(a,b)
special.betaln(a,b) #f_3(a,b)5 勒让德函数
1.连带勒让德函数（多项式）
国内教材对连带勒让德函数大多定义为：
P_{l}^m(x)=(1-x^2)^\frac{m}{2}\frac{\text{d}^m}{\text{d}x^m}P_l(x)\\
而SciPy上的连带勒让德函数定义为：
P_{l}^m(x)=(-1)^m(1-x^2)^\frac{m}{2}\frac{\text{d}^m}{\text{d}x^m}P_l(x)\\
也就是前面多了一个 (-1)^m ，这里需要注意。
SciPy中连带勒让德函数的用法：
special.lpmv(m,l,x)m、l、x的意义与上式相同，m为求导的阶数，l为勒让德多项式的阶数，x为自变量。
2.球谐函数
SciPy中提供的球谐函数是已经归一化了的：
Y_{lm}(\theta,\varphi)=\sqrt{\frac{2l+1}{4\pi}\cdot\frac{(l-|m|)!}{(l+|m|)!}}P_{l}^{|m|}(\cos\theta)e^{im\varphi}\\
special.sph_harm(m,l,varphi,theta)这里需要注意，输入的 m ， n 要保证： |m|\le n ， n\ge0 ，自变量要保证： \varphi\in[0,2\pi] ， \theta\in[0,\pi] 。
3.第一类勒让德函数（多项式）
第一类勒让德函数就是我们常用的勒让德多项式：
P_l(x)=\frac{1}{2^ll!}\frac{\text{d}^l}{\text{d}x^l}(x^2-1)^l\\
SciPy提供的函数为：
special.lpn(l,x)需要注意的是，该函数x的位置不支持输入数组。而且该函数的返回值不太符合直觉。
该函数的返回值是：
([P_0(x), P_1(x), P_2(x), ..., P_l(x)],