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

一維波動方程

對于一根兩端固定、沒有受到任何外力的弦, 若只研究其中的一段, 在不太長的時間 里, 固定端來不及對這段弦產(chǎn)生影響, 則可以認為固定端是不存在的, 弦的長度為無限大。 這種無界 $(?\infty <x<\infty )$ 弦的自由振動由式 $(1)$ 描述。
$$\frac{{?}^{2}u}{?t^2}=a^2\frac{{?}^{2}u}{?x^2}\text{\hspace{1em}(1)}$$
如果保證數(shù)值計算的區(qū)間足夠大, 在一定時間內(nèi), 弦的振動范圍始終沒有超出計算區(qū)間 (或可以近似地這么認為), 那么就能夠放心地使用周期性邊界條件。取 $a=1$, 初始 條件為:
$$u{u}_{t=0}=2\mathrm{sech}(x),{\frac{?u}{?t}|}_{t=0}=0\text{\hspace{1em}(2)}$$
在數(shù)學(xué)物理方法中, 無界弦的自由振動可由行波法求出解析解, 即達朗貝爾公式。 根據(jù)達朗貝爾公式, 從 $t=0$ 開始, $u$ 的初始狀態(tài) $2\mathrm{sech}(x)$ 將分裂為兩個 sech 形的波, 分別向兩邊以速度 $a$ 傳播出去, 即正行波和反行波。下面用傅里葉縉方法求解無界弦 的自由振動問題, 并與達朗貝爾公式的預(yù)測進行比較。首先引入函數(shù) $v$ 對式 $(1)$ 進行降階:
$$\{\begin{array}{l}\frac{?u}{?t}=v\\ \frac{?v}{?t}=a^2\frac{{?}^{2}u}{?x^2}\end{array}\text{\hspace{1em}(3)}$$
對上式等號兩邊做傅里葉變換, 化為偏微分方程組:
$$\{\begin{array}{l}\frac{?\hat{u}}{?t}=\hat{v}\\ \frac{?\hat{v}}{?t}=?a^2{k}^2\hat{u}\end{array}\text{\hspace{1em}(4)}$$

這樣就可以用 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:4subplot(5,2,n)plot(x,usol(p(n),:),'k','LineWidth',1.5),xlabel x,ylabel utitle(['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

計算結(jié)果如圖所示, 初始狀態(tài)的波形分裂成兩半, 并分別向 $x$ 軸正方向和負方向 以速度 $a$ 運動, 這和達朗貝爾公式給出的結(jié)論是一致的。