数学物理方程的求解

• 为了不让我们的文章充斥着难以读懂的乱码，我们还是先上LaTeX比较好

数学物理定解问题的求解方法

• 行波法

• 分离变量法

• 幂级数解

• 格林函数法

• 积分变换

• 包角变换

• 变分法

• 数值解

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import scipy.special as ss
import sympy as sy

特殊函数

• 开始前，一定要仔细检查文献中这些函数的具体定义。文献不同，定义也会小有不同！

Airy Functions

Elliptic Functions

Gamma Functions

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import scipy.special as ss
from scipy.special import factorial
import sympy as sy


fig = plt.figure(figsize=(8,24),tight_layout=True)
fig.suptitle("Gamma")

ax1 = plt.subplot(311)
ax1.set_title("Gamma Function: $\Gamma(z)=\int_0^{\infty}e^{-t}t^{z-1}dt$")

X = np.linspace(0,5,101)
Y = ss.gamma(X)
ax1.plot(X,Y,label="Gamma Function in Real Sets")
plt.legend()


ax2 = plt.subplot(312)
ax2.set_title("Beta Function: $ B(P,Q)=\int_0^1 x^{P-1}(1-x)^{Q-1} dx $")
for p in range(2,6):
    for q in range(2,6):
        Y = ss.beta(p,q)
        ax2.scatter(p,q)
        ax2.text(p+0.1,q-0.1,str(round(Y,3)))
        ax2.set_xlabel("p")
        ax2.set_ylabel("q")
plt.legend()


ax3 = plt.subplot(313)
ax3.set_title("Psi Function $\Psi$")
Y = ss.psi(X)
ax3.plot(X,Y,label="$\Psi function$")
ax3.text(3,-2,"$Euler Const:-\psi(1)=0.577 215 ...$")
plt.pause(0.2)
plt.savefig("0.jpg")

Bessel Funcitons

• v阶贝塞尔方程

• v不为整数

• v任意值

• 汉克尔函数

fig = plt.figure(figsize=(8,8))
fig.suptitle("Bessel")

ax1 = plt.subplot(211)
ax1.set_title("FIRST KIND")

X = np.linspace(0,10,101)
for i in range(5):
    Y = ss.jv(i,X)
    ax1.plot(X,Y,label=str(i)+"th order of first kind Bessel function")
    plt.legend()
ax1.axhline(0)

ax2 = plt.subplot(212)
ax2.set_title("SECOND KIND")
for i in range(5):
    Y = ss.yn(i,X[10:])
    ax2.plot(X[10:],Y,label=str(i)+"th order of second kind Bessel function")
    plt.legend()


plt.pause(0.2)
plt.savefig("0.jpg")

Laguerre Functions

import numpy as np
import scipy.special as ss
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(6,6))
x = np.arange(0,10.01,0.01)
for i in range(5):
    plt.plot(x,ss.laguerre(i)(x),label=str(i))
plt.legend()
plt.savefig("0.jpg")
plt.pause(0.01)

Legendre Functions

还有递推公式：

• 注意：python scipy.special.lpmv 里的连带勒让德函数乘上了一个 (-1)^m

help(ss.legendre)
ss.legendre(6)
poly1d([ 1.44375000e+01,  0.00000000e+00, -1.96875000e+01,  1.60288449e-15,
        6.56250000e+00,  0.00000000e+00, -3.12500000e-01])

help(ss.lpmv)
X = np.linspace(-1,1,102)
Y = ss.lpmv(1,2,X)
plt.plot(X,Y)
f = lambda x:3*x*(1-x**2)**0.5
plt.plot(X,f(X))

以俯仰角为变量的勒让德函数

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import scipy.special as ss
import sympy as sy


fig = plt.figure(figsize=(12,4))

ax1 = plt.subplot(131,projection = "polar")
ax2 = plt.subplot(132,projection = "polar")
ax3 = plt.subplot(133)


def abs_lpmv_polar(m,l,x):
    return abs(((-1)**m)*ss.lpmv(m,l,x))

def lpmv_polar(m,l,x):
    return ((-1)**m)*ss.lpmv(m,l,x)

def lpmv_rect(m,l,x):
    return ((-1)**m)*ss.lpmv(m,l,x)
 
X = np.linspace(0,np.pi,100)

Y1 = abs_lpmv_polar(1,3,np.cos(X))
Y2 = lpmv_polar(1,3,np.cos(X))
Y3 = lpmv_rect(1,3,np.cos(X))

ax1.plot(X,Y1,label="m1l3")
ax1.set_title("abs_lpmv_polar")
ax1.legend()

ax2.plot(X,Y2,label="m1l3")
ax2.set_title("lpmv_polar")
ax2.legend()

ax3.plot(X,Y3,label="m1l3")
ax3.legend()

plt.pause(0.2)

实球谐函数

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import scipy.special as ss
from scipy.special import factorial
import sympy as sy


fig = plt.figure(figsize=(12,4))

ax1 = plt.subplot(131,projection = "3d")
ax2 = plt.subplot(132,projection = "3d")
ax3 = plt.subplot(133,projection = "3d")

def lpmv(m,l,x):
    return ((-1)**m)*ss.lpmv(m,l,x)

def Yfun(l,theta,phi):
    A = ((2*l+1)/4/np.pi)**0.5
    return A * lpmv(0,l,theta) 

Theta = np.linspace(0,np.pi,101)
Phi = np.linspace(0,2*np.pi,201)

Theta,Phi = np.meshgrid(Theta,Phi)

X = Yfun(1,Theta,Phi)*np.sin(Theta)*np.cos(Phi)
Y = Yfun(1,Theta,Phi)*np.sin(Theta)*np.sin(Phi)
Z = Yfun(1,Theta,Phi)*np.cos(Theta)
ax1.plot_surface(X,Y,Z)


X = Yfun(2,Theta,Phi)*np.sin(Theta)*np.cos(Phi)
Y = Yfun(2,Theta,Phi)*np.sin(Theta)*np.sin(Phi)
Z = Yfun(2,Theta,Phi)*np.cos(Theta)
ax2.plot_surface(X,Y,Z)

X = Yfun(3,Theta,Phi)*np.sin(Theta)*np.cos(Phi)
Y = Yfun(3,Theta,Phi)*np.sin(Theta)*np.sin(Phi)
Z = Yfun(3,Theta,Phi)*np.cos(Theta)
ax3.plot_surface(X,Y,Z)

plt.pause(0.2)
plt.savefig("0.jpg")

数学物理方程数值解法

有限元法

多重网格方法

有限差分方法

有限体积法

近似求解的误差估计方法

多尺度计算方法