简介

最开始看到这个名字，我也很激动，终于有个中文姓氏的数学公式了，然鹅司徒卢威是个俄国人，而且司徒卢威完全是音译，就离谱。

司徒卢威函数是下面的非齐次贝赛尔方程的一组解：

$$x^2\frac{{\text{d}}^{2}y}{\text{d}{x}^2}+x\frac{\text{d}y}{\text{d}x}+(x^2-v^2)y=\frac{4(\frac{x}{2}{)}^{v+1}}{\sqrt{\pi }\Gamma (v+\frac{1}{2})}$$

定义为

$$H_v(z)=(z/2{)}^{v+1}\sum _{n=0}^{\infty }\frac{(-1{)}^n(z/2{)}^{2n}}{\Gamma (n+\frac{3}{2})\Gamma (n+v+\frac{3}{2})}$$

图像

当$v=0,1,2,3$时，其函数图像如图所示

绘图函数为

import numpy as np
from scipy.special import struve
import matplotlib.pyplot as plt

x = np.linspace(-10., 10., 1000)
for v in range(4):
    plt.plot(x, struve(v, x), label=f'$H_{i!r}$')

plt.legend(ncol=2)
plt.show()

当$v$取值为非整数时，要求$x$取值必须大于0，可以绘制动图来查看其变化

# 这三个包在后面的程序中不再复述
import matplotlib.animation as animation

fig = plt.figure(figsize=(7,4))
ax = fig.add_subplot(autoscale_on=False, 
    xlim=(0,20),ylim=(0,200))
ax.grid()
plt.tight_layout()
theta_text = ax.text(0.05,0.95,'',
    transform=ax.transAxes)

line, = ax.plot([], [], lw=1)

xs = np.linspace(0., 20., 1000)
def animate(v):
    line.set_data(xs, struve(v, xs))
    theta_text.set_text(f'v={v}')
    return line, theta_text

vs = np.arange(50)/5
ani = animation.FuncAnimation(fig, animate, vs, 
    interval=15, blit=True)
ani.save("struve.gif")
plt.show()

效果为

其他相关函数

scipy中提供了一些与司徒卢威函数联系紧密的函数，包括

func	函数	表达式
modstruve	修正司徒卢威函数	$L_v(x)=-i\exp (-i\pi v/2)H_v(x)$
itstruve0	$H_0$的积分	${\int }_0^x{H}_0(t)\text{d}t$
it2struve0	与$H_0$有关的的积分	${\int }_x^{\infty }\frac{H_0(t)}{t}\text{d}t$
itmodstruve0	修正$H_0$的积分	${\int }_0^x{L}_0(t)\text{d}t$

其中modstruve需要指定$v$，在其取值不同时函数图像为

代码如下

import scipy.special as ss
x = np.linspace(-10., 10., 1000)
for v in range(4):
    plt.plot(x, ss.modstruve(v, x), label=f'$L_{i!r}$')

plt.show()

三个积分函数的图像为

绘图代码如下

funcDct = {
    'integral modified H0':ss.itmodstruve0,
    'integral H0':ss.itstruve0,
    'integral H0(t)/(t)':ss.it2struve0
}

fig = plt.figure()
xs = np.linspace(-10,10,1000)
for i,key in enumerate(funcDct):
    ax = fig.add_subplot(1,3,1+i)
    ys = funcDct[key](xs)
    ax.plot(xs, ys)
    ax.set_title(key)

plt.show()