高效解法：Matlab微分方程谱方法实现

考虑如下的二维波动方程：

可认为上式描述了水的波动过程，u(x,y,t)是水面高度的分布，x、y是空间坐标，t是时间。在|x|=3边界处采用了周期性边界条件，这代表从x=3边界流出的水将从x=-3边界流入，反之亦然。在|y|=1边界处采用了诺依曼边界条件，∂u/∂y|_|y|=1=0代表水不能从|y|=1处流入或流出。

由初始条件可知∂u/∂y|_{|y|=1,t=0≈0}，因此可将边界条件∂u/∂y|_|y|=1=0等价地写为∂(∂u/∂y)/∂t|_|y|=1=0。又由它的连续性，所以可交换求导顺序，等价写为∂(∂u/∂t)/∂y|_|y|=1=0。然后，引入函数v(x,y,t)=∂u(x,y,t)/∂t，将式（5-74）中的∂²/∂t²降为∂/∂t，得到：

用ode45计算上式中的∂/∂t。此外，由于|x|=3处存在周期性边界条件，所以可使用谱求导矩阵计算x方向上的导数。对于|y|=1边界处的诺依曼边界条件，这里使用切比雪夫求导矩阵计算y方向上的导数，并根据式（5-69）处理∂v/∂y|_|y|=1=0。代码如下：

程序5-20

主程序代码如下：(www.xing528.com)

文件wave_tank.m代码如下：

计算结果如图5-24所示，边界|y|=1像“水池壁”一样将“水”挡住。方程式（5-74）来源于文献[2]，这里使用了比文献[2]更精确的方法，有兴趣的读者可以比较一下两种方法的差异。

图5-24 二维波动方程的解，|x|=3边界处采用了周期性边界条件，|y|=1边界 处采用了诺依曼边界条件∂u/∂y|_|y|=1=0