谱方法学习笔记-上(超详细)

谱方法学习笔记📒

谱方法学习笔记-下(超详细)
声明：鉴于CSDN使用 $KaTeX$ 渲染公式，$\text{KATEX}$ 与 $LaTeX$ 不同，不支持直接的交叉引用命令，如\label和\eqref。$\text{KATEX}$专注于数学公式的呈现，而不提供完整的文档生成功能。因此，为了保证显示正常，本文中所有公式均已删除编号和交叉引用。由此可能导致一些可读性上的不便，请读者谅解。

傅立叶谱方法

对于在 $(-\infty ,\infty )$ 有定义且绝对可积、并在任一有限区间上满足狄利克莱条件的函数 $u(x)$, 傅里叶变换（Fourier transform）及其逆变换（inverse Fourier transform） 定义为:
$$\hat{u}(k)={\int }_{-\infty }^{\infty }u(x){\mathrm{e}}^{-\mathrm{i}kx}\mathrm{d}x$$

$$u(x)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\hat{u}(k){\mathrm{e}}^{\mathrm{i}kx}\mathrm{d}k$$

上述定义式给出了一个傅里叶变换对 (Fourier transform pair), 本专栏中将它们记为 $\hat{u}(k)=F[u(x)]$ 和 $u(x)=F^{-1}[\hat{u}(k)]$ .实际上, 傅里叶变换对的定义并不是唯一的, 两个定义式中的系数可以随意修改, 只要它们的积为 $1/2\pi$ 即可.此外, 定义式中的 ${\mathrm{e}}^{-\mathrm{i}kx}$ 和 ${\mathrm{e}}^{\mathrm{i}kx}$ (称为积分变换的核) 也可以互换.容易知道, $u(x)$ 先后经过傅里叶变换及其逆变换仍将得到它本身, 即: u ( x ) = F − 1 { F [ u ( x ) ] } u(x)=F^{-1}\{F[u(x)]\} u(x)=F−1{F[u(x)]} .

通常来讲, 若不对 “傅里叶变换” 加任何限定, 那么它指的就是连续傅里叶变换, 也 就是针对定义在无限区间内的连续函数 $u(x)$ 的傅里叶变换, 这是在理想条件下的数学定义. 在实际应用中, 尤其是计算机的信号采样、信号处理当中, 信号是离散的、有限的, 离散傅里叶变换（discrete Fourier transform) 就是针对这一情况提出的.在 Matlab 中, 对于序列 $u_1,\dots ,u_j,\dots ,u_N$ 的离散傅里叶变换及逆变换定义为:
$${\hat{u}}_k=\sum _{j=1}^N{u}_j{\mathrm{e}}^{\frac{-2\pi (j-1)(k-1)\mathrm{i}}{N}},k=1,\dots ,N$$

$$u_j=\frac{1}{N}\sum _{k=1}^N{\hat{u}}_k{\mathrm{e}}^{\frac{2\pi (j-1)(k-1)\mathrm{i}}{N}},j=1,\dots ,N$$

同样, 上述定义中的归一化系数也可以有其他选择, 但它们的乘积必须为 $1/N$ .如果将序列 $u_1,\dots ,u_j,\dots ,u_N$ 看作等间隔空间（时间）点上的信号幅度值, 那么经过离散傅里叶变换得到的序列 ${\hat{u}}_1,\dots ,{\hat{u}}_k,\dots ,{\hat{u}}_N$ 就是其相应的频谱信息.通过定义式 % \eqref{3-4} 和 % \eqref{3-5} 容易得到 ${\hat{u}}_k={\hat{u}}_{k+N}$ 和 $u_j=u_{j+N}$, 这就是说, 离散傅里叶变换已经隐含了周期性边界条件.

在 Matlab 中，fft 和 ifft 函数实现基于快速傅里叶变换算法的一维离散傅里叶变换及逆变换。fft 函数的一般调用语法为：

fft(u)

若u为一向量，则 fft 返回其离散傅里叶变换的结果。若 u为一矩阵，则 fft 返回矩阵中每一列的离散傅里叶变换。若要计算矩阵 u 中每一行的离散傅里叶变换，只需调用 fft(u,[],2)即可，它速度比 fft(u.').'更快一些。ifft 函数的调用语法与 fft 函数完全一样。

若输入 fft 函数的空域信号序列 $u_1,...,u_j,...,u_N$ 在 $x$ 轴上的坐标间隔是 $L/N$ , 相应的横坐标 $x$ 可认为是：

$$-L/2,-L/2+L/N,\dots ,L/2-2L/N,L/2-L/N$$

注意，fft 函数输出的序列 ${\hat{u}}_1,...,{\hat{u}}_k,...,{\hat{u}}_N$ 所对应的频率并不是按照从小到大的顺序排列的，而是非负频率对应于前半部分序列，负频率对应于后半部分序列，即：

$$0,2\pi /L,4\pi /L,\dots ,(N-4)\pi /L,(N-2)\pi /L,-N\pi /L,-(N-2)\pi /L,\dots ,-4\pi /L,-2\pi /L$$

因此 Matlab 提供了 fftshift 函数，用于调整 fft 输出序列的顺序。也就是说，fftshift(fft(u)) 所对应的频率才是按照递增顺序排列的：

$$-N\pi /L,-(N-2)\pi /L,\dots ,-2\pi /L,0,2\pi /L,\dots ,(N-4)\pi /L,(N-2)\pi /L$$

空域坐标与频域坐标的关系如下表所示，不难发现，空域上的两点间隔和频域区间长度成反比，乘积为 $2\pi$, 频域上的两点间隔和空域区间长度的关系亦如此。

	两点间隔	区间长度
空域坐标	$L/N$	$L$
频域坐标	$2\pi /L$	$2\pi N/L$

此外还有 ifftshift 函数用于进行与 fftshift 函数相反的操作。实际上，fftshift 函数的作用仅是让 fft 的输出结果变得更易于观察，若还需要用 ifft 函数将频谱信息还原回空域信号，就必须在使用 fftshift 之后再用 ifftshift 函数将序列的顺序恢复，否则得不到正确结果。如果不是为了分析频谱信息，而是用下文即将介绍的傅里叶谱方法处理导数和微分问题，则大可不必使用 fftshift 和 ifftshift 函数调整序列的顺序，以减少代码冗余、 节约时间。

基于快速傅里叶变换算法的二维离散傅里叶变换及逆变换由 fft2、ifft2 函数实现，一般调用语法为：

fft2(u)

它返回矩阵 u 的二维傅里叶变换，等效于 fft(fft(u),[],2)。ifft2 与 fft2 的调用语法一样，就不再赘述。

fftshift 函数的调用语法为：

fftshift(u)

若u为一向量，它将 u 左右两半互换并返回。若 u 为一矩阵，则它将 u 的第 1、3 象限互换以及第 2、4 象限互换并返回。此外，若要互换矩阵 u 的上下两半，可以调用 fftshift(u,1)。要互换矩阵 u 的左右两半，只需调用 fftshift(u,2)。ifftshift 函数与 fftshift 函数功能完全相反， 调用语法相同，这里不再重复。

求导与傅里叶谱方法

对于 $F\left[u^{\prime }(x)\right]$, 由傅里叶变换的定义和分部积分法, 得到:
$$F\left[u^{\prime }(x)\right]={\int }_{-\infty }^{\infty }{u}^{\prime }(x){\mathrm{e}}^{-\mathrm{i}kx}\mathrm{d}x={u(x){\mathrm{e}}^{-\mathrm{i}kx}|}_{-\infty }^{\infty }-{\int }_{-\infty }^{\infty }u(x)(-\mathrm{i}k){\mathrm{e}}^{-\mathrm{i}kx}\mathrm{d}x$$
当 $\mid x\mid \to \infty$ 时, $u(x)\to 0$, 则:
$$F\left[u^{\prime }(x)\right]=\mathrm{i}k{\int }_{-\infty }^{\infty }u(x){\mathrm{e}}^{-\mathrm{i}kx}\mathrm{d}x=\mathrm{i}kF[u(x)]$$
类似地, 可以得到:
$$F\left[u^{(n)}(x)\right]=(\mathrm{i}k{)}^{n}F[u(x)]$$
其中 $u^{(n)}(x)$ 代表 $u(x)$ 的 $n$ 阶导数.上式的意义在于：函数的求导运算在傅里叶变换的作用下, 可转化为相对简单的代数运算, 即: u ( n ) ( x ) = F − 1 { ( i k ) n F [ u ( x ) ] } u^{(n)}(x)=F^{-1}\left\{(\mathrm{i} k)^n F[u(x)]\right\} u(n)(x)=F−1{(ik)nF[u(x)]} .正是基于此原理, 傅里叶谱方法 (Fourier spectral method) 利用傅里叶变换将偏微分方程中空域 (space domain）或时域（time domain）上的求导运算简化为频域 ( spectral domain) 上的代数运算, 求解后再通过傅里叶逆变换得到空域或时域上的结果.在代码层面上, Matlab 提供的快速傅里叶变换函数 fft、逆变换函数 ifft 以及强大的矩阵运算能力也为简洁、优雅地实现傅里叶谱方法奠定了基础.

如前所述，fft 输出结果在频域上的顺序是颠倒的，经过 fft 函数变换后的序列在 $k$ 轴上所对应的坐标是：
$$0,2\pi /L,4\pi /L,...,(N-4)\pi /L,(N-2)\pi /L,-N\pi /L,-(N-2)\pi /L,...,-4\pi /L,-2\pi /L$$
利用傅里叶谱方法求导时，在代码中通常使用顺序颠倒的频域坐标，以求序列 u 的一阶导数数值解为例：ifft(i*k.*fft(u))。其中，变量 k 是由 Matlab 语句 (2*pi/L)*[0:N/2-1 -N/2:-1] 生成的顺序颠倒的频域坐标，这样做是为了免去调用 fftshift 和 ifftshift 的步骤。

下面用实例说明傅里叶谱方法求 $u(x)={\mathrm{e}}^{\sin (\pi x)}$ 的 $1、2、3$ 阶导数的过程, 并与解析解

$$u^{\prime }(x)=\pi \cos (\pi x){\mathrm{e}}^{\sin (\pi x)}$$
$$u^{\prime \prime }(x)={\pi }^2{\mathrm{e}}^{\sin (\pi x)}\left[{\cos }^2(\pi x)-\sin (\pi x)\right]$$
$$u^{\prime \prime \prime }(x)={\pi }^3{\mathrm{e}}^{\sin (\pi x)}\cos (\pi x)\left[{\cos }^2(\pi x)-3\sin (\pi x)-1\right]$$

相比较, 得到的最大误差分别在 $10^{-14}、10^{-13}、10^{-11}$ 数量级, 如图所示.
代码如下:

clear all; close all;
L=2; N=32;
x=L/N*[-N/2:N/2-1];
k=2*pi/L*[0:N/2-1 -N/2:-1];
u=exp(sin(pi*x)); ut=fft(u);
%导数的精确解
du_exact(:,1)=pi*cos(pi*x).*u;
du_exact(:,2)=pi^2*(cos(pi*x).^2-sin(pi*x)).*u;
du_exact(:,3)=pi^3*cos(pi*x).*(cos(pi*x).^2-3*sin(pi*x)-1).*u;
%谱方法求导
du_Fourier(:,1)=ifft(i*k.*ut);
du_Fourier(:,2)=ifft((i*k).^2.*ut);
du_Fourier(:,3)=ifft((i*k).^3.*ut);
error=max(abs(du_exact-du_Fourier));
%画图
labels={'u''(x)','u''''(x)','u''''''(x)'};
for n=1:4
    subplot(2,2,n)
    if n==1
        plot(x,u,'k','LineWidth',1.5)
        xlabel x, ylabel u(x), title('u(x)=e^{sin(\pix)}')
    else
        plot(x,du_exact(:,n-1),'k',x,du_Fourier(:,n-1),'.r' ...
            ,'MarkerSize',13,'LineWidth',1.5)
        title(['Error_{max}=' num2str(error(n-1))])
        xlabel x, ylabel(labels(n-1))
    end
end

傅立叶谱方法的步聚

待求解的偏微分方程的普遍形式为:
$$\frac{{\partial }^{n}u}{\partial t^n}=Lu+N(u)$$
其中, $u(x,t)$ 为 $x、t$ 的函数, $L$ 代表线性算符 (linear operator), $N(u)$ 为非线性项 (nonlinear terms ).已知初始条件为 $u\left(x,t_0\right)$, 在周期性边界条件下求某一时刻的 $u(x,t)$ .对于这样的问题, 总是可以通过变量代换的办法将等号左边对 $t$ 的 $n$ 阶导数降为 1 阶, 所以下面以等号左边是 $\partial u/\partial t$ 的偏微分方程为例进行讨论, 这与式 % \eqref{3-11} 降阶后的偏微分方程组求解方法大同小异.
在 $x$ 域上对以下方程做傅里叶变换:
$$\frac{\partial u}{\partial t}=Lu+N(u)$$
得到:
$$\frac{\partial \hat{u}}{\partial t}=\alpha (k)\hat{u}+F[N(u)]$$
其中, $\hat{u}(k,t)$ 代表 $u(x,t)$ 在 $x$ 域上的傅里叶变换.为了展现对线性算符 $L$ 及非线性项 $N(u)$ 做傅里叶变换的细节, 这里以

$$L=a\cdot {\partial }^2/\partial x^2+b\cdot \partial /\partial x+c$$

$$N(u)=u^3+u^3\cdot {\partial }^{2}u/\partial x^2+f(x)\cdot \partial u/\partial x$$

为例进行分析, $a、b$ 和 $c$ 分别为常数.

对线性部分, 有 $F[Lu]=\left[a(\mathrm{i}k{)}^2+b(\mathrm{i}k)+c\right]\hat{u},$ 则式 % \eqref{3-13} 中的 $\alpha (k)=a(\mathrm{i}k{)}^2+b(\mathrm{i}k)+c=$ $-a{k}^2+\mathrm{i}bk+c$ .

对非线性部分, 需要将 $u、\partial u/\partial x$ 和 ${\partial }^{2}u/\partial x^2$ 写为 $F^{-1}[\hat{u}]、F^{-1}[\mathrm{i}k\hat{u}]$ 和 $F^{-1}\left[(\mathrm{i}k{)}^2\hat{u}\right]$, 则有
$$\begin{array}{rl}F[N(u)]& =F\left[u^3+u^3\cdot {\partial }^{2}u/\partial x^2+f(x)\cdot \partial u/\partial x\right]\\ & =F\left\{F^{-1}[\hat{u}{]}^3+F^{-1}[\hat{u}{]}^3\cdot F^{-1}\left[(\mathrm{i}k{)}^2\hat{u}\right]+f(x)\cdot F^{-1}[\mathrm{i}k\hat{u}]\right\}\end{array}$$
这样, 式 % \eqref{3-12} 中 $u(x,t)$ 对 $x$ 的偏导数就都通过傅里叶变换及逆变换简化为式 % \eqref{3-13} 中 $\hat{u}(k,t)$ 和 $k$ 的代数运算了, 然后再将 $\hat{u}$ 和 $k$ 离散化, 偏微分方程就简化成了常微分方程组.数值计算 $\partial \hat{u}(k,t)/\partial t$ 最直接的方法就是调用 Matlab 的 ode 系列函数, 优先选择 ode 45 函数, 最后用 ifft 函数将频域上的计算结果 $\hat{u}(k,t)$ 变换回待求的 $u(x,t)$ .

综上, 利用傅里叶谱方法数值计算偏微分方程 (组) 步骤如图（一维情况下的傅里叶谱方法计算过程）所示, 具体过程如下:

• 对于形如式 % \eqref{3-11} 的偏微分方程 (组), 通过变量代换将其中的 ${\partial }^{n}u/\partial t^n$ 降为 1 阶 导数, 若 $n=1$, 略过此步.

• 在 $x$ 域上对偏微分方程 (组) 做傅里叶变换, 则线性项中的 ${\partial }^{n}u/\partial x^n$ 直接变为 $(\mathrm{i}k{)}^n\hat{u}$, 并利用 $u=F^{-1}[\hat{u}]、{\partial }^{n}u/\partial x^n=F^{-1}\left[(\mathrm{i}k{)}^n\hat{u}\right]$ 代换非线性项中的 $u、{\partial }^{n}u/\partial x^n$, 得到关于 $\hat{u}(k,t)$ 的微分方程 (组).

• 用时间步进法 (如龙格-库塔法等) 数值计算离散化的关于 $\hat{u}(k,t)$ 的微分方程组, 默认边界条件为周期性边界条件, 初始条件 $\hat{u}\left(k,t_0\right)=F\left[u\left(x,t_0\right)\right]$ .

• 将上一步得到的结果从频域变换回空域, 即: $u(x,t)=F^{-1}[\hat{u}(k,t)]$ .离散傅里叶变换及逆变换由 fft、ifft 函数实现.

傅里叶谱方法求解基本偏微分方程—一维波动方程

对于一根两端固定、没有受到任何外力的弦, 若只研究其中的一段, 在不太长的时间里, 固定端来不及对这段弦产生影响, 则可以认为固定端是不存在的, 弦的长度为无限大。 这种无界 $(-\infty <x<\infty )$ 弦的自由振动由式 $(1)$ 描述。
$$\frac{{\partial }^{2}u}{\partial t^2}=a^2\frac{{\partial }^{2}u}{\partial x^2}$$
如果保证数值计算的区间足够大, 在一定时间内, 弦的振动范围始终没有超出计算区间 (或可以近似地这么认为), 那么就能够放心地使用周期性边界条件。取 $a=1$, 初始条件为:
$${u|}_{t=0}=2\mathrm{sech}(x),{\frac{\partial u}{\partial t}|}_{t=0}=0$$
在数学物理方法中, 无界弦的自由振动可由行波法求出解析解, 即达朗贝尔公式。 根据达朗贝尔公式, 从 $t=0$ 开始, 长度为 $L=80$ 的 $u$ 的初始状态 $2\mathrm{sech}(x)$ 将分裂为两个 sech 形的波, 分别向两边以速度 $a=1$ 经过 $T=20$ 秒传播出去, 即正行波和反行波。下面用傅里叶谱方法求解无界弦的自由振动问题, 并与达朗贝尔公式的预测进行比较。首先引入函数 $v$ 对式 % \eqref{111} 进行降阶:
$$\{\begin{array}{l}\frac{\partial u}{\partial t}=v\\ \frac{\partial v}{\partial t}=a^2\frac{{\partial }^{2}u}{\partial x^2}\end{array}$$
对上式等号两边做傅里叶变换（连续）, 化为偏微分方程组:
$$\{\begin{array}{l}\frac{\partial \hat{u}}{\partial t}=\hat{v}\\ \frac{\partial \hat{v}}{\partial t}=-a^2{k}^2\hat{u}\end{array}$$

注意，以上操作均在连续系统下进行，未做离散近似。

接下来开始构建数值格式，首先，对区间 $[-\frac{L}{2},\frac{L}{2}]$ 做空间剖分，得到空域信号序列 $u_1,...,u_j,...,u_N$ 在 $x$ 轴上的坐标间隔是 $L/N$ , 相应的横坐标 $x$ 可认为是：
$$-L/2,-L/2+L/N,\dots ,L/2-2L/N,L/2-L/N$$

此时可使用 Matlab 中 fft 函数做离散傅立叶变换，需要注意的是，fft 函数输出的序列 ${\hat{u}}_1,...,{\hat{u}}_k,...,{\hat{u}}_N$ 所对应的频率并不是按照从小到大的顺序排列的，而是非负频率对应于前半部分序列，负频率对应于后半部分序列，即相应的横坐标 $k$ 可认为是：
$$0,2\pi /L,4\pi /L,\dots ,(N-4)\pi /L,(N-2)\pi /L,-N\pi /L,-(N-2)\pi /L,\dots ,-4\pi /L,-2\pi /L$$

因此 Matlab 提供了 fftshift 函数，用于调整 fft 输出序列的顺序。也就是说，fftshift(fft(u)) 所对应的频率才是按照递增顺序排列的：
$$-N\pi /L,-(N-2)\pi /L,\dots ,-2\pi /L,0,2\pi /L,\dots ,(N-4)\pi /L,(N-2)\pi /L$$
此外还有 ifftshift 函数用于进行与 fftshift 函数相反的操作。实际上，fftshift 函数的作用仅是让 fft 的输出结果变得更易于观察，若还需要用 ifft 函数将频谱信息还原回空域信号，就必须在使用 fftshift 之后再用 ifftshift 函数将序列的顺序恢复，否则得不到正确结果。如果不是为了分析频谱信息，而此处目的是通过傅里叶谱方法处理微分问题，则大可不必使用 fftshift 和 ifftshift 函数调整序列的顺序，而是在代码中通常使用顺序颠倒的频域坐标，以减少代码冗余、 节约时间。

以求序列 u 的二阶导数数值解为例：ifft((i*k).^2.*fft(u))。其中，变量 k 是由 Matlab 语句 (2*pi/L)*[0:N/2-1 -N/2:-1] 生成的顺序颠倒的频域坐标，这样做是为了免去调用 fftshift 和 ifftshift 的步骤。

这样就可以用 ode45 求解了, 详细代码如下:

主程序代码如下:

clear all; close all;

L=80;N=256;
x=L/N*[-N/2:N/2-1];
k=(2*pi/L)*[0:N/2-1 -N/2:-1].';
% 初始条件
u=2*sech(x);ut=fft(u);
vt=zeros(1,N);uvt=[ut vt];
% 求解
a=1;t=0:0.5:20;
[t,uvtsol]=ode45('wave1D',t,uvt,[],N,k,a);
usol=ifft(uvtsol(:,1:N),[],2);
% 画图
p=[1 11 21 41];
for n=1:4
    subplot(5,2,n)
    plot(x,usol(p(n),:),'k','LineWidth',1.5),xlabel x,ylabel u
    title(['t=' num2str(t(p(n)))]),axis([-L/2 L/2 0 2])
end
subplot(5,2,5:10)
waterfall(x,t,usol),view(10,45)
xlabel x,ylabel t,zlabel u,axis([-L/2 L/2 0 t(end) 0 2])

文件 wave1D.m 代码如下:

function duvt=wave1D(t,uvt,dummy,N,k,a)
ut=uvt(1:N);vt=uvt(N+[1:N]);
duvt=[vt;-a^2*(k).^2.*ut];
end

计算结果如图所示, 初始状态的波形分裂成两半, 并分别向 $x$ 轴正方向和负方向 以速度 $a$ 运动, 这和达朗贝尔公式给出的结论是一致的。

傅里叶谱方法求解基本偏微分方程—二维波动方程

将一维波动方程中的一维无界弦自由振动方程推广到二维空间上, 就得到了描述无界 $(-\infty <x,y<\infty )$ 弹性薄膜的波动方程:
$$\frac{{\partial }^{2}u}{\partial t^2}=a^2\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u$$
取 $a=1$, 初始条件为:
$${u|}_{t=0}={\mathrm{e}}^{-20\left[(x-0.4{)}^2+(y+0.4{)}^2\right]}+{\mathrm{e}}^{-20\left[(x+0.4{)}^2+(y-0.4{)}^2\right]},{\frac{\partial u}{\partial t}|}_{t=0}=0$$
可以这样理解上述初始条件的物理意义: 两手抓住弹性薄膜的两个位置, 分别提起, 使薄膜上形成两个峰, 在 $t=0$ 时刻突然松手。根据生活常识可以预料到, 这两个位置的薄 膜将来回振动, 与此同时, 产生的波向四周传播, 而且波与波会在相遇处叠加。
为便于求解, 引入函数 $v$ 对式 % \eqref{222} 进行降阶, 得:
$$\{\begin{array}{l}\frac{\partial u}{\partial t}=v\\ \frac{\partial v}{\partial t}=a^2\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u\end{array}$$

对上式等号两边做傅里叶变换, 得到常微分方程组:
$$\{\begin{array}{l}\frac{\partial \hat{\tilde{u}}}{\partial t}=\hat{\hat{v}}\\ \frac{\partial \hat{\hat{v}}}{\partial t}=-a^2\left(k_x^2+k_y^2\right)\hat{\hat{u}}\end{array}$$
接下来用 ode45 求解即可, 代码如下:

主程序代码如下:

clear all; close all;

L=4;N=64;
x=L/N*[-N/2:N/2-1];y=x;
kx=(2*pi/L)*[0:N/2-1 -N/2:-1];ky=kx;
[X,Y]=meshgrid(x,y);
[kX,kY]=meshgrid(kx,ky);
K2=kX.^2+kY.^2;
% 初始条件
u=exp(-20*((X-0.4).^2+(Y+0.4).^2))+exp(-20*((X+0.4).^2+(Y-0.4).^2));
ut=fft2(u);vt=zeros(N);uvt=[ut(:); vt(:)];
% 求解
a=1;t=[0 0.25 0.5 1];
[t,uvtsol]=ode45('wave2D',t,uvt,[],N,K2(:),a);
% 画图
for n=1:4
    subplot(2,2,n)
    mesh(x,y,ifft2(reshape(uvtsol(n,1:N^2),N,N))),view(10,45)
    title(['t=' num2str(t(n))]),axis([-L/2 L/2 -L/2 L/2 0 1])
    xlabel x,ylabel y,xlabel x,zlabel u
end

文件 wave2D.m 代码如下:

function duvt=wave2D(t,uvt,dummy,N,K2,a)
ut=uvt(1:N^2);vt=uvt(N^2+[1:N^2]);
duvt=[vt;-a^2*K2.*ut];
end

程序输出结果如图所示, 它反映了弹性薄膜上的波向四周传播的过程。