傅里葉譜方法求解基本偏微分方程—二維波動方程

二維波動方程

將一維波動方程中的一維無界弦自由振動方程推廣到二維空間上, 就得到了描述無界 $(?\infty <x,y<\infty )$ 彈性薄膜的波動方程:
$$\frac{{?}^{2}u}{?t^2}=a^2\left(\frac{{?}^2}{?x^2}+\frac{{?}^2}{?y^2}\right)u\text{\hspace{1em}(1)}$$
取 $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{?u}{?t}|}_{t=0}=0\text{\hspace{1em}(2)}$$
可以這樣理解上述初始條件的物理意義: 兩手抓住彈性薄膜的兩個位置, 分別提起, 使薄膜上形成兩個峰, 在 $t=0$ 時刻突然松手。根據(jù)生活常識可以預(yù)料到, 這兩個位置的薄 膜將來回振動, 與此同時, 產(chǎn)生的波向四周傳播, 而且波與波會在相遇處疊加。
為便于求解, 引入函數(shù) $v$ 對式 $(1)$ 進(jìn)行降階, 得:
$$\{\begin{array}{l}\frac{?u}{?t}=v\\ \frac{?v}{?t}=a^2\left(\frac{{?}^2}{?x^2}+\frac{{?}^2}{?y^2}\right)u\end{array}\text{\hspace{1em}(3)}$$

對上式等號兩邊做傅里葉變換, 得到常微分方程組:
$$\{\begin{array}{l}\frac{?\hat{\tilde{u}}}{?t}=\hat{\hat{v}}\\ \frac{?\hat{\hat{v}}}{?t}=?a^2\left(k_x^2+k_y^2\right)\hat{\hat{u}}\end{array}\text{\hspace{1em}(4)}$$
接下來用 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:4subplot(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

文件 wave1D.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

程序輸出結(jié)果如圖所示, 它反映了彈性薄膜上的波向四周傳播的過程。