Python科学计算：偏微分方程1

首先，我们来看初边值问题：伯格斯方程：

假设函数是定义在上的函数，且满足：

右侧第一项表示自对流，第二项则表示扩散，在许多物理过程中，这两种效应占据着主导地位，为了固定一个特定的解，我们对其施加一个初始条件：

以及一个或者多个边值条件：

由上面的三个式子所组成的问题被称为初边值问题(IBVP)，如果我们同时设置a为-inf，b为 inf，那么我们会得到一个初值问题(IVP)

这里主要介绍两个比较常用的方法：

1、有限差分空间导数：我们选用有限组等距值来表示区间

其中步长，那么我们根据导数的定义，可以得到：

以及：

这样的话，我们就可以用前向欧拉方法来计算，但是，这样的话，我们就要进行稳定性分析(之前的博客里面是有的)，为了保证显式时间步进是稳定的，就要确保：,其中C是一个同阶常数，那么从这里我们实际上可以看得到：显式差分的时间步长收到了空间步长的制约，这是显式差分的一个弊病。

2、周期问题的谱技术空间导数方法：说实话，我之前一直很想用伪谱法，但是因为本科期间没有学过复数和傅里叶变换等，所以一直不敢动手做，今天有这个机会，那就得好好学一学。

谱方法是有限差分法的一个有用的替代方法，仍然使用伯格斯方程进行讲解：首先考虑空间域被映射到频率域上的一个特殊情况，假设的x具有周期，如果用来表示

相对于x的傅里叶变换，因为：的傅里叶变换是，那么通过傅里叶变换，可以恢复的空间导数，所以，对于伯格斯方程，需要一个傅里叶变换和两个傅里叶逆变换来恢复等式右侧的两个导数，需要把其变成一个光谱算法：

假设在区间的n个等距点上用函数值表示每一个t时刻的，那么就可以构造离散傅里叶变换(DFT)，其广义上是的傅里叶级数展开的前N项，为每个项乘以适当的乘数，然后计算逆dft。

如果我们知道是平滑的，也就是说，存在任意多个x导数并且有界，那么对于任意大的k，DFT的截断误差是，如果，就大致可以返回阶的误差，有限差分法的则需要N~

来达到同样的精度，如果N只有小的素数因子，比如说，那么DFT以及其逆变换可以采用快速傅里叶变换FFT来计算，需要次的运算，关键是，既是周期性的，也是平滑性的，如果

，那么周期延拓将会是不连续的，并且吉布斯现象将破坏所有这些估计的精确性。

scipy中有一个非常有用的函数diff，如果输入的是一个numpy数组，表示在上均匀间隔的值，则函数返回一个与u形状相同的数组，包含相同x值的order阶导数的值。

看一个书上给的例子：计算导数值

import numpy as np
from scipy.fftpack import diffdef fd(u):"返回2*dx有限差分"ud=np.empty_like(u)ud[1:-1]=u[2:]-u[:-2]ud[0]=u[1]-u[-1]ud[-1]=u[0]-u[-2]return ud
for N in [4,8,16,32,64,128,256]:dx=2.0*np.pi/Nx=np.linspace(0,2.0*np.pi,N,endpoint= False)u=np.exp(np.sin(x))du_ex=np.cos(x)*udu_sp=diff(u)du_fd=fd(u)/(2.0*dx)err_sp=np.max(np.abs(du_sp-du_ex))err_fg=np.max(np.abs(du_fd-du_ex))print("N=%d,err_sp=%.4e,err_sp=%.4e"% (N,err_fg,err_sp))

图1

从图一可以看出，点数每增加一倍，有限差分法的误差的无限范数（最大绝对值）就减少大约4倍，光谱误差岁每个加倍呈现平方的增大，直到N非常大，N>30，这种快速的误差减小方式通常被称作指数收敛。

我们再来看一个空间周期问题的IVP：

考虑线性对流：

其精确解是

图2

结果图形说明了伪谱法的优点：即使是较粗糙的网格间距也能够产生平滑的结果。

import numpy as np
from scipy.fftpack import diff
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#6-19建立问题
def u_exact(t,x):"Exact solution"return np.exp(np.sin(x-2*np.pi*t))
def rhs(u,t):"return rhs."return -2.0*np.pi*diff(u)
N=32
x=np.linspace(0,2*np.pi,N,endpoint=False)
u0=u_exact(0,x)
t_initial=0.0
t_final=64*np.pi
t=np.linspace(t_initial,t_final,101)
#20行用来求解
sol=odeint(rhs,u0,t,mxstep=5000)
#可视化
fig=plt.figure()
ax=fig.add_subplot(1,1,1, projection='3d')
t_gr,x_gr=np.meshgrid(x,t)
ax.plot_surface(t_gr,x_gr,sol,alpha=0.5)
ax.elve,ax.azim=47,-137
ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_zlabel("u")
plt.show()
u_ex=u_exact(t[-1],x)
err=np.abs(np.max(sol[-1,:]-u_ex))
print("with %d fourier nodes the final error = %g"%(N,err))

再，再看看，理解理解。