譜方法學(xué)習(xí)筆記📒

譜方法學(xué)習(xí)筆記-上(超詳細(xì))
聲明：鑒于CSDN使用 $KaTeX$ 渲染公式，$\text{KATE?X}$ 與 $LaTeX$ 不同，不支持直接的交叉引用命令，如\label和\eqref。$\text{KATE?X}$專注于數(shù)學(xué)公式的呈現(xiàn)，而不提供完整的文檔生成功能。因此，為了保證顯示正常，本文中所有公式均已刪除編號(hào)和交叉引用。由此可能導(dǎo)致一些可讀性上的不便，請(qǐng)讀者諒解。

傅里葉譜方法求解基本偏微分方程—二維 Schnakenberg 模型

Schnakenberg模型是一種描述化學(xué)反應(yīng)動(dòng)力學(xué)的數(shù)學(xué)模型，旨在研究化學(xué)反應(yīng)中自組織現(xiàn)象和空間模式的形成。

該模型最早由德國(guó)化學(xué)家Heinrich Otto Wieland提出，后由德國(guó)數(shù)學(xué)家Theodor Schnakenberg在其博士論文中進(jìn)行了推導(dǎo)和解析，并得到了該模型的名字。

Schnakenberg模型的意義在于，它可以模擬一系列反應(yīng)中出現(xiàn)的自組織現(xiàn)象，例如化學(xué)波、化學(xué)斑圖、化學(xué)分叉等。這些現(xiàn)象是由于反應(yīng)物濃度分布在空間上的不均勻性導(dǎo)致的，通過(guò)Schnakenberg模型可以更好地理解這些現(xiàn)象，并為實(shí)驗(yàn)提供一定的指導(dǎo)和驗(yàn)證。

此外，Schnakenberg模型還可以應(yīng)用于其他領(lǐng)域，例如生物學(xué)、物理學(xué)、生態(tài)學(xué)等。在這些領(lǐng)域中，該模型可以用于研究許多自組織現(xiàn)象，例如細(xì)胞分裂、病毒傳播、生態(tài)系統(tǒng)演化等。

這個(gè)模型描述的是兩種化學(xué)物質(zhì)之間的反應(yīng)過(guò)程。具體來(lái)說(shuō)，化學(xué)物質(zhì) $u$ 會(huì)被消耗掉，同時(shí)生成化學(xué)物質(zhì)$v$。反應(yīng)速率的大小取決于$u$ 和$v$的濃度，而擴(kuò)散系數(shù)則描述了物質(zhì)在空間中的擴(kuò)散速度。

Schnakenberg模型有許多有趣的特性和行為。例如，在某些參數(shù)范圍內(nèi)，這個(gè)模型的解可以顯示出空間上的分布和時(shí)間上的演化，形成有序的斑圖（patterns）。這些斑圖包括靜態(tài)和動(dòng)態(tài)的形態(tài)，它們的形狀和數(shù)量取決于模型參數(shù)的取值。這種有序結(jié)構(gòu)的出現(xiàn)和演化，可以通過(guò)Schnakenberg模型來(lái)解釋和理解實(shí)驗(yàn)室中觀察到的許多化學(xué)反應(yīng)和生物學(xué)現(xiàn)象。具體來(lái)說(shuō)，當(dāng)模型的參數(shù)取值適當(dāng)時(shí)，Schnakenberg模型的解會(huì)形成一個(gè)包含不同濃度級(jí)別的周期性結(jié)構(gòu)。這種周期性結(jié)構(gòu)被稱為Turing結(jié)構(gòu)，是由Alan Turing在1952年提出的一種理論模型。Turing結(jié)構(gòu)是許多自組織系統(tǒng)中的一種普遍現(xiàn)象，例如斑馬條紋、蝴蝶翅膀上的斑點(diǎn)等。
在Schnakenberg模型中，Turing結(jié)構(gòu)的出現(xiàn)是由反應(yīng)和擴(kuò)散之間的相互作用所導(dǎo)致的。具體來(lái)說(shuō)，當(dāng)擴(kuò)散系數(shù)足夠大時(shí)，化學(xué)物質(zhì)可以在空間上擴(kuò)散，使得濃度分布均勻。但當(dāng)反應(yīng)速率較慢時(shí)，這些化學(xué)物質(zhì)在局部區(qū)域內(nèi)聚集起來(lái)，并且在一定的空間尺度上發(fā)生反應(yīng)。這個(gè)過(guò)程會(huì)導(dǎo)致局部區(qū)域內(nèi)的濃度發(fā)生變化，從而形成周期性的結(jié)構(gòu)。
需要注意的是，Turing結(jié)構(gòu)的出現(xiàn)是非線性動(dòng)力學(xué)的一個(gè)重要特征。這種結(jié)構(gòu)的產(chǎn)生和演化需要考慮多個(gè)因素之間的相互作用，因此通常需要使用數(shù)值計(jì)算和數(shù)學(xué)分析的方法來(lái)研究其性質(zhì)和行為。

需要注意的是，Schnakenberg模型是一種簡(jiǎn)化的模型，它假設(shè)了化學(xué)物質(zhì)的濃度在空間和時(shí)間上均勻分布，沒(méi)有考慮一些具體的生物或化學(xué)體系的細(xì)節(jié)和復(fù)雜性。因此，這個(gè)模型只能作為一種基本的參考模型，而不是具體體系的精確描述。

近年來(lái)，Schnakenberg模型被廣泛應(yīng)用于不同領(lǐng)域的研究中，包括流體力學(xué)、生物學(xué)、材料科學(xué)等。這個(gè)模型的理論和數(shù)值分析方法也在不斷地發(fā)展和完善，為我們深入理解和探究復(fù)雜動(dòng)態(tài)系統(tǒng)的性質(zhì)和行為提供了有力的工具和平臺(tái)。

斑圖（pattern）是一類普遍存在于自然界、在時(shí)間或空間上具有某種規(guī)律的非均勻宏觀結(jié)構(gòu)。反應(yīng)-擴(kuò)散系統(tǒng) (reaction-diffusion system) 是斑圖理論中研究得最為廣 泛的系統(tǒng), 它起源于化學(xué)反應(yīng)系統(tǒng), 但又不局限于此, 還廣泛應(yīng)用于生物學(xué)、物理學(xué)、醫(yī)學(xué)、金融學(xué)等。Schnakenberg 模型是反應(yīng)-擴(kuò)散系統(tǒng)中的一個(gè)有趣的模型, 數(shù)學(xué)形式如下:
$$?\begin{array}{l}\frac{?u}{?t}=\left(\frac{{?}^2}{?x^2}+\frac{{?}^2}{?y^2}\right)u+\gamma \left(a?u+u^{2}v\right)\\ \frac{?v}{?t}=d\left(\frac{{?}^2}{?x^2}+\frac{{?}^2}{?y^2}\right)v+\gamma \left(b?u^{2}v\right)\end{array}$$
其中, $u$ 和 $v$ 可看做兩種化學(xué)反應(yīng)物質(zhì)的濃度, $x、y$ 為空間坐標(biāo), $t$ 為時(shí)間, 、 $\gamma$ 為常數(shù)。對(duì)式 (3-38) 做二維傅里葉變換, 可將其轉(zhuǎn)化為偏微分方程組:

$$\{\begin{array}{l}\frac{?\hat{u}}{?t}=?\left(k_x^2+k_y^2\right)\hat{\hat{u}}+\gamma ?F\left\{a?F^{?1}[\hat{u}]+F^{?1}[\hat{u}{]}^2{F}^{?1}[\hat{v}]\right\}\\ \frac{?\hat{v}}{?t}=?d\left(k_x^2+k_y^2\right)\hat{v}+\gamma ?F\left\{b?F^{?1}[\hat{u}{]}^2{F}^{?1}[\hat{v}]\right\}\end{array}$$

參數(shù)取值為: $a=0.1,b=0.8,d=26,\gamma =100$ 。為了得到靶型波, 將初始條件設(shè)置為: $u$ 在 $x?y$ 平面原點(diǎn)處為 1 , 在其他位置為 $0,v$ 在整個(gè) $x?y$ 平面上均為 1 。利用傅里葉譜方法 求解該模型的代碼如下。

主程序代碼

clear all; close all;
L=16; N=64;
kx=2*pi/L*[0:N/2-1 -N/2:-1]; ky=kx;
[kX,kY]=meshgrid(kx,ky);
K2=kX.^2+kY.^2;
%初始條件
u=zeros(N); u(N/2,N/2)=1; v=ones(N);
ut=fft2(u); vt=fft2(v);
uvt=[ut(:); vt(:)];
%求解
a=0.1; b=0.8; d=26; gamma=100; t=[0:0.1:0.3];
[t,uvtsol]=ode45('schnakenberg',t,uvt,[],K2,N,gamma,a,b,d);
%畫圖
for n=1:4subplot(2,2,n)gca=pcolor(ifft2(reshape(uvtsol(n,1:N^2),N,N))); axis offset(gca,'LineStyle','none'), shading interptitle(['t=' num2str(t(n))]), axis('square'),
%     colormap('gray')
end

程序輸出的結(jié)果如圖所示，這與化學(xué)實(shí)驗(yàn)中觀察到的靶型波一致。

Matlab 將計(jì)算結(jié)果制作成 gif 動(dòng)畫

引言

求解包含時(shí)間的偏微分方程 (組) 將得到隨著時(shí)間變化的數(shù)值結(jié)果, 把這樣的數(shù)據(jù)制作成 gif 動(dòng)畫并結(jié)合到幻燈片中, 在畢業(yè)答辯、小組討論、課堂教學(xué)等場(chǎng)合有著廣泛的應(yīng)用。生動(dòng)的彩色 gif 動(dòng)畫具有很強(qiáng)的表現(xiàn)力, 令人刮目相看, 大大提高了報(bào)告人所講述理論結(jié)果的直觀性、生動(dòng)性、觀賞性。

函數(shù)介紹

生成 gif 動(dòng)畫主要用到 $4$ 個(gè)函數(shù): getframe、frame2im、rgb2ind、imwrite。

1. getframe 函數(shù)的一般調(diào)用形式為: F=getframe(h), 其作用是截取句柄為 h 的窗口內(nèi)的一幀圖像。

2. frame2im 函數(shù)的作用是把一幀截圖轉(zhuǎn)為圖像數(shù)據(jù)。

3. rgb2ind 函數(shù)的作用是將 RGB 圖像轉(zhuǎn)換為索引圖像, 一般調(diào)用形式為: [X, map]= rgb2ind(RGB,n) 。其中, X、map 分別為轉(zhuǎn)換后的圖像數(shù)據(jù)和顏色表數(shù)據(jù), RGB 為轉(zhuǎn)換前的圖像數(shù)據(jù), n 指定 map 中的顏色數(shù)。

4. imwrite 函數(shù)的作用是將圖像數(shù)據(jù)寫入圖像文件, 一般調(diào)用形式為: imwrite(X,map,filename,fmt,Param1,Val1,Param2,Val2...)。其中, X、map 意義同上, filename 為文件名, fmt 為文件格式,Param1,Val1,Param2,Val2... 為若干可選參數(shù)及其取值。如：參數(shù) LoopCount 為動(dòng)畫的循環(huán)播放次數(shù), 這里設(shè)為 inf, 即無(wú)窮大。參數(shù) DelayTime 為每幀間隔時(shí)間, 單位秒。參數(shù) WriteMode 為寫入文件的模式, 有覆蓋 overwrite (默認(rèn)) 和追加 append 兩種選擇。

Matlab源代碼

生成 gif 動(dòng)畫的示例代碼如下:

---

clear
clc
close
x=-1:0.02:1;y=x;
[X,Y]=meshgrid(x,y);
filename='test.gif';
for a=1:10u=a*exp(-10*(X.^2+Y.^2));mesh(x,y,u),axis([-1 1 -1 1 0 10]),drawnowim=frame2im(getframe(gcf));[A,map]=rgb2ind(im,256);if a==1%先以覆蓋模式寫入指定的gif文件imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.1);else%再以追加模式將每一幀寫入gif文件imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.1);end
end

---

運(yùn)行代碼之后在當(dāng)前目錄下生成 gif 文件, 該動(dòng)畫顯示了一個(gè)三維高斯函數(shù)的峰值逐漸增大的過(guò)程。

為了更加直觀的感受求解 Schnakenberg 模型得到的靶型波，此處結(jié)合上述操作將計(jì)算結(jié)果制作成如下 gif 動(dòng)畫。

修改后的程序源代碼：

clear all; close all;
L=16; N=64;
kx=2*pi/L*[0:N/2-1 -N/2:-1]; ky=kx;
[kX,kY]=meshgrid(kx,ky);
K2=kX.^2+kY.^2;
%初始條件
u=zeros(N); u(N/2,N/2)=1; v=ones(N);
ut=fft2(u); vt=fft2(v);
uvt=[ut(:); vt(:)];
%求解
a=0.1; b=0.8; d=26; gamma=100; t=[0:0.01:1];
[t,uvtsol]=ode45('schnakenberg',t,uvt,[],K2,N,gamma,a,b,d);
%畫圖
filename = 'output.gif'; % 設(shè)置保存的gif文件名
for n = 1:length(t)gca=pcolor(ifft2(reshape(uvtsol(n,1:N^2),N,N))); axis offset(gca,'LineStyle','none'), shading interptitle(['t=' num2str(t(n))]), axis('square'),% colormap('gray')% 捕捉當(dāng)前子圖并轉(zhuǎn)換為gif圖像幀frame = getframe(gcf);im=frame2im(getframe(gcf));[A,map]=rgb2ind(im,256);% 寫入gif圖像文件if n == 1%先以覆蓋模式寫入指定的gif文件imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.01);else%再以追加模式將每一幀寫入gif文件imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.01);end
end

譜求導(dǎo)矩陣

譜求導(dǎo)矩陣的導(dǎo)出和應(yīng)用

在生產(chǎn)實(shí)踐過(guò)程中, 某些函數(shù)關(guān)系的具體表達(dá)式是未知的, 只能通過(guò)實(shí)驗(yàn)測(cè)量得到該函數(shù)上的一些離散的數(shù)據(jù)點(diǎn).所以, 人們希望通過(guò)這些離散的數(shù)據(jù)點(diǎn)構(gòu)造一個(gè)已知函數(shù)來(lái)近似地描述并替代未知函數(shù).此外, 盡管有些函數(shù)的表達(dá)式是已知的, 但由于太過(guò)復(fù)雜繁瑣, 不便對(duì)其進(jìn)行理論計(jì)算和數(shù)值分析, 所以也有必要構(gòu)造一個(gè)簡(jiǎn)單函數(shù)來(lái)替代它.

上述問(wèn)題用數(shù)學(xué)語(yǔ)言表達(dá)為: 已知在區(qū)間 $[a,b]$ 上的 $N$ 個(gè)位置 $x_1,x_2,\dots ,x_N$ 的函數(shù)值 $u_1$, $u_2,\dots ,u_N$, 求一個(gè)函數(shù) $p(x)$ 通過(guò)所有已知點(diǎn) $\left(x_1,u_1\right),\left(x_2,u_2\right),\dots ,\left(x_N,u_N\right)$, 即:

$$u_j=p\left(x_j\right),j=1,2,\dots ,N$$

其中, 用 $p(x)$ 近似未知函數(shù)的方法稱為插值法 (interpolation), $p(x)$ 稱為插值函數(shù).下面介紹用譜方法確定插值函數(shù) $p(x)$ 的過(guò)程, 并將其用于計(jì)算離散數(shù)據(jù)的各階導(dǎo)數(shù)、偏微分方程 (組) 的數(shù)值求解.

譜方法插值

考慮在區(qū)間 $[0,2\pi ]$ 上、具有周期性邊界條件的插值問(wèn)題.在此區(qū)間上間距為 $h$ 的 $N$ 個(gè)位置 $x_1,x_2,\dots ,x_N$ 對(duì)應(yīng)的函數(shù)值為 $u_1,u_2,\dots ,u_N$ .其中, $x_j=jh,h=2\pi /N$, 即:

$$\frac{\pi }{h}=\frac{N}{2}$$

在推導(dǎo)前先做 3 點(diǎn)說(shuō)明:

(1) 選擇區(qū)間 $[0,2\pi ]$ 是為了方便討論, 在該區(qū)間上得出的結(jié)論同樣適用于對(duì)其平移得到的其他區(qū)間 (如 $[?\pi ,\pi ]$ ), 并可以通過(guò)乘以縮放因子方便地將結(jié)論轉(zhuǎn)化到其他長(zhǎng)度非 $2\pi$的任意區(qū)間上 (如 $[?L,L])$ .

(2) 所謂周期性邊界條件, 是指 $N$ 個(gè)函數(shù)值 $u_1,u_2,\dots ,u_N$ 可等效地看做是無(wú)窮多個(gè)函數(shù)值 $\dots ,u_{?1},u_0,u_1,\dots ,u_{N?1},u_N,u_{N+1},\dots$ 的一部分, 并存在 $u_{j+mN}=u_j$ 的關(guān)系 ( $m$ 為任意整數(shù)).這里把討論的函數(shù)當(dāng)做周期為 $2\pi$ 的周期函數(shù)處理.

(3) $N$ 的奇偶會(huì)導(dǎo)致接下來(lái)的推導(dǎo)細(xì)節(jié)有所差異, 但過(guò)程是相似的, 所以此處只分析 $N$ 為偶數(shù)的情況.

若對(duì)序列 $u_1,u_2,\dots ,u_N$ 做離散傅里葉變換, 那么空域（時(shí)域）上的間隔 $h$ 決定了頻域上的區(qū)間為 $[?\pi /h,\pi /h]$ .由式 % \eqref{4-2} , 該區(qū)間也可寫為 $[?N/2,N/2]$ .本部分中的離散傅里葉變換對(duì)定義為:

$$\begin{array}{r}{\hat{u}}_k=h\sum _{j=1}^N{\mathrm{e}}^{?\mathrm{i}k{x}_j}{u}_j,k=?\frac{N}{2}+1,\dots ,\frac{N}{2}{u}_j=\frac{1}{2\pi }\sum _{k=?N/2+1}^{N/2}{\mathrm{e}}^{\mathrm{i}k{x}_j}{\hat{u}}_k,j=1,\dots ,N\end{array}\text{\hspace{1em}?}$$

這與前面給出的 Matlab 中離散傅里葉變換對(duì)的定義略有不同, 但實(shí)質(zhì)是一樣的.注意到式 % \eqref{4-4} 中的頻率分量并不是完全對(duì)稱的, $k$ 中有 ,卻沒(méi)有 $?N/2$, 這在求解插值函數(shù)時(shí)會(huì)引起一些小小的問(wèn)題.所以, 令 ${\hat{u}}_{?N/2}={\hat{u}}_{N/2}$, 并重新定義離散傅里葉逆變換為:

$$u_j=\frac{1}{2\pi }\sum _{k=?N/2}^{N/2}{}^{\prime }{\mathrm{e}}^{\mathrm{i}k{x}_j}{\hat{u}}_k,j=1,\dots ,N$$

其中, “ ${\Sigma }^{\prime }$ 代表求和時(shí)在 $k=\pm N/2$ 的項(xiàng)上乘以 $1/2$ .需要強(qiáng)調(diào)的是, 式 % \eqref{4-3} 和式 % \eqref{4-4} 仍然是離散傅里葉變換對(duì)的定義, 式 % \eqref{4-5} 僅是用于確定插值函數(shù) $p(x)$ 的.求 $p(x)$ 時(shí)需要把式 % \eqref{4-5} 中的 $x_j=jh$ 推廣到 $[0,2\pi ]$ 上的任意實(shí)數(shù) $x$, 即:

$$p(x)=\frac{1}{2\pi }\sum _{k=?N/2}^{N/2}{}^{\prime }{\mathrm{e}}^{\mathrm{i}kx}{\hat{u}}_k,x\in [0,2\pi ]$$

確定插值函數(shù)的步驟是這樣的：先用式 % \eqref{4-3} 將序列 $u_1,u_2,\dots ,u_N$ 變換為 ${\hat{u}}_{?N/2+1}$, ${\hat{u}}_{?N/2+2},\dots ,{\hat{u}}_{N/2}$, 令 ${\hat{u}}_{?N/2}={\hat{u}}_{N/2}$, 再根據(jù)序列 ${\hat{u}}_{?N/2},{\hat{u}}_{?N/2+1},\dots ,{\hat{u}}_{N/2}$ 通過(guò)式 % \eqref{4-6} 得到 $p(x)$ .這樣得到的插值函數(shù) $p(x)$ 可用來(lái)求序列 $u_1,u_2,\dots ,u_N$ 在 $x_j$ 處的各階導(dǎo)數(shù), 即:

$${\frac{{?}^{n}p(x)}{?x^n}|}_{x=x_j}$$

上面介紹的是通過(guò)譜方法計(jì)算插值函數(shù) $p(x)$ 來(lái)估算序列 $u_1,u_2,\dots ,u_N$ 在 $x_j$ 處的導(dǎo)數(shù)的基本原理.接下來(lái), 為了把上述過(guò)程轉(zhuǎn)化為方便的矩陣運(yùn)算, 采用如下思路: 首先求出周期 $\delta$ 函數(shù)的插值函數(shù) $S_N(x)$, 然后將任意序列 $u_1,u_2,\dots ,u_N$ 寫為周期 $\delta$ 函數(shù)的線性組合, 進(jìn)而可把它的插值函數(shù) $p(x)$ 寫為 $S_N(x)$ 的線性組合, 最后找到 $p(x)$ 和 $S_N(x)$ 在 $x_1,x_2,\dots ,x_N$ 處導(dǎo)數(shù)的關(guān)系并寫為矩陣形式, 給出針對(duì)任意序列 $u_1,u_2,\dots ,u_N$ 的譜求導(dǎo)矩陣.

譜求導(dǎo)矩陣

周期 $\delta$ 函數(shù)定義為:

$${\delta }_j=\{\begin{array}{ll}1,& (j\%N=0)\\ 0,& (j\%N\ne 0)\end{array}$$

其中, 百分號(hào) “ $\%$ ” 代表求余運(yùn)算, 周期 $\delta$ 函數(shù)在 $j=mN$ ( $m$ 為任意整數(shù)) 時(shí)取值為 1 ,其他情況為 0 , 它所對(duì)應(yīng)的橫坐標(biāo)為 $x_j=jh$ .利用式 % \eqref{4-3} 求周期 $\delta$ 函數(shù)的離散傅里葉變換, 結(jié)果為一常數(shù) $h$ :

$${\hat{\delta }}_k=h\sum _{j=1}^N{\mathrm{e}}^{?\mathrm{i}k{x}_j}{\delta }_j=h,k=?\frac{N}{2}+1,\dots ,\frac{N}{2}$$

再利用式 % \eqref{4-6} 求周期 $\delta$ 函數(shù)的插值函數(shù):

$$\begin{array}{rl}p(x)& =\frac{h}{2\pi }\sum _{k=?N/2}^{N/2}{}^{\prime }{\mathrm{e}}^{\mathrm{i}kx}\\ & =\frac{h}{2\pi }\left(\frac{1}{2}\sum _{k=?N/2}^{N/2?1}{\mathrm{e}}^{\mathrm{i}kx}+\frac{1}{2}\sum _{k=?N/2+1}^{N/2}{\mathrm{e}}^{\mathrm{i}kx}\right)\\ & =\frac{h}{2\pi }\left(\frac{1}{2}{\mathrm{e}}^{?\frac{\mathrm{i}}{2}x}\sum _{k=?N/2+1/2}^{N/2?1/2}{\mathrm{e}}^{\mathrm{i}kx}+\frac{1}{2}{\mathrm{e}}^{\frac{\mathrm{i}}{2}x}\sum _{k=?N/2+1/2}^{N/2?1/2}{\mathrm{e}}^{\mathrm{i}kx}\right)\\ & =\frac{h}{2\pi }\cos (x/2)\sum _{k=?N/2+1/2}^{N/2?1/2}{\mathrm{e}}^{\mathrm{i}kx}\\ & =\frac{h}{2\pi }\cos (x/2)\frac{{\mathrm{e}}^{\mathrm{i}(?N/2+1/2)x}?{\mathrm{e}}^{\mathrm{i}(N/2+1/2)x}}{1?{\mathrm{e}}^{\mathrm{i}x}}\\ & =\frac{h}{2\pi }\cos (x/2)\frac{{\mathrm{e}}^{?\mathrm{i}(N/2)x}?{\mathrm{e}}^{\mathrm{i}(N/2)x}}{{\mathrm{e}}^{?\mathrm{i}x/2}?{\mathrm{e}}^{\mathrm{i}x/2}}\\ & =\frac{h}{2\pi }\cos (x/2)\frac{\sin (Nx/2)}{\sin (x/2)}\end{array}$$

由式 % \eqref{4-2} , 最終得到的周期 $\delta$ 函數(shù)的插值函數(shù)為周期 sinc 函數(shù) $S_N$ :

$$S_N(x)=\frac{\sin (\pi x/h)}{(2\pi /h)\tan (x/2)}$$

可以證明, ${S_N(x)|}_{x\to 0}=1$ .下圖（曲線為周期 $\mathrm{sinc}$ 函數(shù), 黑點(diǎn)為周期 $\delta$ 函數(shù), 上下兩圖分別對(duì)應(yīng) $N=4$ 和 $N=16$） 顯示了周期 $\delta$ 函數(shù)以及它的插值函數(shù)一一周期 sinc 函數(shù) $S_N(x)$, 二者的周期均為 $2\pi$, 無(wú)論 $N$ 取值為 4 還是 16 , 后者都精確、巧妙地經(jīng)過(guò)了前者的所有離散點(diǎn).

若將區(qū)間 $[0,2\pi ]$ 上的序列 $u_1,u_2,\dots ,u_N$ 寫為周期 $\delta$ 函數(shù)的線性疊加, 即:
$$u_j=\sum _{m=1}^N{u}_m{\delta }_{j?m}$$

那么, 將其中的周期 $\delta$ 函數(shù)替換為周期 $\mathrm{sinc}$ 函數(shù) $S_N$, 就得到序列 $u_1,u_2,\dots ,u_N$ 的插值函數(shù) $p(x)$ :

$$p(x)=\sum _{m=1}^N{u}_m{S}_N\left(x?x_m\right)$$

注意到 $p(x)$ 就是 $S_N(x)$ 的線性組合,如果用 $p(x)$ 來(lái)估算序列 $u_1,u_2,\dots ,u_N$ 在 $x_j=jh$ 處的 1 階導(dǎo)數(shù), 則有:

$$p^{\prime }\left(x_j\right)=\sum _{m=1}^N{u}_m{S}_N^{\prime }\left(x_j?x_m\right)$$

其中, $x_j?x_m=(j?m)h$ .將上式寫為矩陣形式:

$$\left(\begin{array}{c}{p}^{\prime }\left(x_1\right)\\ p^{\prime }\left(x_2\right)\\ p^{\prime }\left(x_3\right)\\ ?\\ p^{\prime }\left(x_N\right)\end{array}\right)=\left(\begin{array}{ccccc}{S}_N^{\prime }(0)& S_N^{\prime }(?h)& S_N^{\prime }(?2h)& ?& S_N^{\prime }((1?N)h)\\ S_N^{\prime }(h)& S_N^{\prime }(0)& S_N^{\prime }(?h)& & \\ S_N^{\prime }(2h)& S_N^{\prime }(h)& S_N^{\prime }(0)& & \\ ?& & & ?& \\ S_N^{\prime }((N?1)h)& & & & \end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ ?\\ u_N\end{array}\right)$$

因此, 只需在向量 ${\left(u_1,u_2,\dots ,u_N\right)}^{\mathrm{T}}$ 上乘以一個(gè) $N$ 階方陣即可得到由它的譜方法插值函數(shù) $p(x)$的導(dǎo)數(shù)組成的向量 ${\left(p^{\prime }\left(x_1\right),p^{\prime }\left(x_2\right),\dots ,p^{\prime }\left(x_N\right)\right)}^{\mathrm{T}}$, 這個(gè) $N$ 階方陣就是譜求導(dǎo)矩陣 (spectral differentiation matrix).下文用 ${\mathbf{D}}_N$ 來(lái)表示 1 階 $N\times N$ 譜求導(dǎo)矩陣, ${\mathbf{D}}_N^{(n)}$ 表示 $n$ 階 $N\times N$ 譜求導(dǎo)矩陣.周期 $\mathrm{sinc}$ 函數(shù) $S_N(x)$ 在 $x_j=jh$ 處的 1 階導(dǎo)數(shù)為:

$$S_N^{\prime }\left(x_j\right)=\{\begin{array}{cc}0,& (j\%N=0)\\ (?1{)}^j/2?\cot (jh/2),& (j\%N\ne 0)\end{array}$$

將其代入到式 % \eqref{4-15} 中的 $N$ 階方陣, 得到 1 階譜求導(dǎo)矩陣:

$${\mathbf{D}}_N=\left(\begin{array}{ccccc}0& & & & ?\frac{\cot (1h/2)}{2}\\ ?\frac{\cot (1h/2)}{2}& ?& & ?& \frac{\cot (2h/2)}{2}\\ \frac{\cot (2h/2)}{2}& ?& & ?\frac{\cot (3h/2)}{2}& \\ ?\frac{\cot (3h/2)}{2}& & ?& & ?\\ ?& ?& ?& \frac{\cot (1h/2)}{2}& \\ \frac{\cot (1h/2)}{2}& & & & 0\end{array}\right)$$

類似地, 還可以得到高階譜求導(dǎo)矩陣.周期 $\mathrm{sinc}$ 函數(shù) $S_N(x)$ 在 $x_j=jh$ 處的 2 階導(dǎo)數(shù)為:

$$S_N^{\prime \prime }\left(x_j\right)=\{\begin{array}{cc}?\frac{{\pi }^2}{3h^2}?\frac{1}{6},& (j\%N=0)\\ ?\frac{(?1{)}^j{\csc }^2(jh/2)}{2},& (j\%N\ne 0)\end{array}$$

那么, 2 階譜求導(dǎo)矩陣為:

$$D_N^{(2)}=\left(\begin{array}{ccc}?& ?& \\ ?& ?\frac{1}{2}{\csc }^2\left(\frac{2h}{2}\right)& \\ ?& \frac{1}{2}{\csc }^2\left(\frac{1h}{2}\right)& \\ & ?\frac{{\pi }^2}{3h^2}?\frac{1}{6}& \\ & \frac{1}{2}{\csc }^2\left(\frac{1h}{2}\right)& ?\\ & ?\frac{1}{2}{\csc }^2\left(\frac{2h}{2}\right)& ?\\ & ?& ?\end{array}\right)$$

周期 sinc 函數(shù) $S_N(x)$ 在 $x_j=jh$ 處的 3 階導(dǎo)數(shù)為:

$$S_N^m\left(x_j\right)=\{\begin{array}{cc}0,& (j\%N=0)\\ (?1{)}^j\cot \left(\frac{jh}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{jh}{2}\right)?\frac{{\pi }^2}{2h^2}\right],& (j\%N\ne 0)\end{array}$$

同樣得到 3 階譜求導(dǎo)矩陣:
$$D_N^{(3)}=(\begin{array}{cccc}0& & & ?\cot \left(\frac{1h}{2}\right)\left[\frac{3}{4}{\cos }^2\left(\frac{1h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]\\ ?\cot \left(\frac{1h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{1h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]& ?& ?& \cot \left(\frac{2h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{2h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]\\ \cot \left(\frac{2h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{2h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]& & ?& ?\cot \left(\frac{3h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{3h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]\\ ?\cot \left(\frac{3h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{3h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]& ?& & ?\\ ?& ?& ?& \cot \left(\frac{1h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{1h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]\\ \cot \left(\frac{1h}{2}\right)\left[\frac{3}{4}{\csc }^2\left(\frac{1h}{2}\right)?\frac{{\pi }^2}{2h^2}\right]& & & 0\end{array})$$

構(gòu)造 $n$ 階譜求導(dǎo)矩陣 ${\mathbf{D}}_N^{(n)}$ 的一般方法為:

(1) 求周期 sinc 函數(shù) $S_{N}(x)$ 在 $x_{j}=j h$ 處的 $n$ 階導(dǎo)數(shù) $S_{N}^{(n)}\left(x_{j}\right)$ .(2) $\boldsymbol{D}_{N}^{(n)}$ 的第 1 列為 $\left(S_{N}^{(n)}\left(x_{N}\right), S_{N}^{(n)}\left(x_{1}\right), S_{N}^{(n)}\left(x_{2}\right), \ldots, S_{N}^{(n)}\left(x_{N-1}\right)\right)^{\mathrm{T}}$, 第 2 列為 $\left(S_{N}^{(n)}\left(x_{N-1}\right)\right.$, $\left.S_{N}^{(n)}\left(x_{N}\right), S_{N}^{(n)}\left(x_{1}\right), \ldots, S_{N}^{(n)}\left(x_{N-2}\right)\right)^{\mathrm{T}}$, 第 3 列為 $\left(S_{N}^{(n)}\left(x_{N-2}\right), S_{N}^{(n)}\left(x_{N-1}\right), S_{N}^{(n)}\left(x_{N}\right), \ldots, S_{N}^{(n)}\left(x_{N-3}\right)\right)^{\mathrm{T}} \cdots \cdots$依此類推.

此外, 周期函數(shù) $S_N^{(n)}(x)$ 的奇偶性是由 $n$ 的奇偶決定的, 這在構(gòu)造 ${\mathbf{D}}_N^{(n)}$ 時(shí)可以產(chǎn)生一些便利:

$${\frac{{?}^n{S}_N(x)}{?x^n}|}_{x=x_j}={(?1{)}^n\frac{{?}^n{S}_N(x)}{?x^n}|}_{x=x_{N?j}}$$

當(dāng) $n$ 為奇數(shù)時(shí), 必有:

$${\frac{{?}^n{S}_N(x)}{?x^n}|}_{x=x_j}=0,j\%N=0$$

當(dāng) $n$ 為偶數(shù)時(shí), $S_N^{(n)}(x)$ 在 $x=x_j$ 處 $(j\%N=0)$ 的取值是無(wú)窮小/無(wú)窮小的形式, 需要用洛必達(dá)法則求它的極限, 即:

$${\frac{{?}^n{S}_N(x)}{?x^n}|}_{x=x_j}=\underset{x\to 0}{\lim }\frac{{?}^n{S}_N(x)}{?x^n},j\%N=0$$

人工計(jì)算 $S_N(x)$ 的 $n$ 階導(dǎo)數(shù)的工作量隨著階數(shù) $n$ 的增大而顯著增加, 比較省時(shí)省力的方法是利用 Matlab 的符號(hào)運(yùn)算功能, 調(diào)用 diff 函數(shù)和 limit 函數(shù)求導(dǎo)、求極限, 并結(jié)合 toeplitz 函數(shù), 可實(shí)現(xiàn)程序自動(dòng)生成 ${\mathbf{D}}_N^{(n)}$, 有興趣的讀者可以自行嘗試.

由于譜求導(dǎo)矩陣 ${\mathbf{D}}_N^{(n)}$ 是在計(jì)算區(qū)間長(zhǎng)度為 $2\pi$ 的前提下構(gòu)造的, 所以若實(shí)際的計(jì)算區(qū)間長(zhǎng)度為 $L$, 則必須在 ${\mathbf{D}}_N^{(n)}$ 上乘以縮放系數(shù)才可保證結(jié)果正確, 即: $(2\pi /L{)}^n{\mathbf{D}}_N^{(n)}$ .

下面給出使用譜求導(dǎo)矩陣 對(duì) $u(x)={\mathrm{e}}^{\sin (\pi x)}$ 求 階導(dǎo)數(shù), 并與精確解

$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]$

比較的實(shí)例.代碼如下:

clear all; close all;
L=2; N=32;
x=L/N*[-N/2:N/2-1]'; 
%構(gòu)造譜求導(dǎo)矩陣
h=2*pi/N; column=[0 0.5*(-1).^(1:N-1).*cot((1:N-1)*h/2)]';
D=(2*pi/L)*toeplitz(column,column([1 N:-1:2]));
column=[-pi^2/(3*h^2)-1/6 -0.5*(-1).^(1:N-1)./sin(h*(1:N-1)/2).^2];
D2=(2*pi/L)^2*toeplitz(column);
column=[0 (-1).^(1:N-1).*cot((1:N-1)*h/2).*( ...-pi^2/(2*h^2)+3/4*csc((1:N-1)*h/2).^2)]';
D3=(2*pi/L)^3*toeplitz(column,column([1 N:-1:2]));
%導(dǎo)數(shù)的精確解
u=exp(sin(pi*x));
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;
%譜方法求導(dǎo)
du_Fourier(:,1)=D*u;
du_Fourier(:,2)=D2*u;
du_Fourier(:,3)=D3*u;
error=max(abs(du_exact-du_Fourier));
%畫圖
labels={'u''(x)','u''''(x)','u''''''(x)'};
for n=1:4subplot(2,2,n)if n==1plot(x,u,'k','LineWidth',1.5)xlabel x, ylabel u(x), title('u(x)=e^{sin(\pix)}')elseplot(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

如下圖（用譜求導(dǎo)矩陣計(jì)算 $u(x)={\mathrm{e}}^{\sin (\pi x)}$ 的 階導(dǎo)數(shù)（點(diǎn)）并與精確解（曲線）比較）所示, 曲線為精確的 和 $u^{\prime \prime \prime }(x)$, 點(diǎn)為利用譜求導(dǎo)矩陣計(jì)算得到的 和 $u^{\prime \prime \prime }\left(x_j\right)$ .在 $N$ 的取值僅為 32 的情況下, 用譜求導(dǎo)矩陣計(jì)算 $u(x)={\mathrm{e}}^{\sin (\pi x)}$ 的 階導(dǎo)數(shù)與精確解的最大誤差分別在 數(shù)量級(jí).

用譜求導(dǎo)矩陣求解偏微分方程的步驟

與之前一樣, 待求解的偏微分方程的普遍形式為:

$$\frac{{?}^{n}u}{?t^n}=Lu+N(u)$$

其中, $u(x,t)$ 為 $x、t$ 的函數(shù), $L$ 代表線性算符, $N(u)$ 為非線性項(xiàng).通過(guò)函數(shù)代換可將等號(hào)左邊的 $n$ 階導(dǎo)數(shù) ${?}^n/?t^n$ 降到 1 階, 所以下面分析式 % \eqref{4-26} 的求解過(guò)程, 這與式 % \eqref{4-25} 降階后得到的方程組的解法是一樣的.

$$\frac{?u}{?t}=Lu+N(u)$$

用譜求導(dǎo)矩陣數(shù)值求解式 % \eqref{4-26} 的步驟如下:

(1) 在 $x$ 軸上的計(jì)算區(qū)間內(nèi), 將 $x$ 離散化為 $N$ 個(gè)等間距的位置 $\mathbf{x}={\left(x_1,x_2,\dots ,x_N\right)}^{\mathrm{T}}$, 相應(yīng)地, 將函數(shù) $u$ 離散化為 $N$ 維向量 $\mathbf{u}={\left(u_1,u_2,\dots ,u_N\right)}^{\mathrm{T}}$ .

(2) 在等號(hào)右邊,將線性算符 $L$ 和非線性項(xiàng) $N(u)$ 里所有對(duì)函數(shù) $u$ 的求導(dǎo)運(yùn)算替換為譜求導(dǎo)矩陣與向量 $\mathbf{u}$ 的乘運(yùn)算, 即: ${?}^{n}u/?x^n\to {\mathbf{D}}_N^{(n)}\mathbf{u}$, 并將 $L$ 和 $N(u)$ 都寫成矩陣形式, 得到形如式 % \eqref{4-27} 的微分方程組.注意, 若計(jì)算區(qū)間長(zhǎng)度不是 $2\pi$, 則必須在 ${\mathbf{D}}_N^{(n)}$ 上乘以縮放系數(shù).

(3) 用時(shí)間步進(jìn)法（歐拉法、龍格-庫(kù)塔法等）數(shù)值計(jì)算式 % \eqref{4-27} 等號(hào)左邊的 $?/?t$,在周期性邊界條件下,得到不同 $t$ 處的向量 $\mathbf{u}$ .

$$\frac{?}{?t}\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ ?\\ u_N\end{array}\right)=\left(\begin{array}{ccccc}{L}_{11}& L_{12}& L_{13}& ?& L_{1N}\\ L_{21}& L_{22}& L_{23}& ?& L_{2N}\\ L_{31}& L_{32}& L_{33}& ?& L_{3N}\\ ?& ?& ?& ?& ?\\ L_{N1}& L_{N2}& L_{N3}& ?& L_{NN}\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ ?\\ u_N\end{array}\right)+N\left[\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ ?\\ u_N\end{array}\right)\right]$$

針對(duì)步驟 2 , 這里以 $L=a?{?}^2/?x^2+b??/?x+c、N(u)=u^3+u^3?{?}^{2}u/?x^2+f(x)??u/?x$ 為例進(jìn)行說(shuō)明, 和 $c$ 分別為常數(shù).對(duì)線性算符 $L$, 除了將 ${?}^n/?x^n$ 替換為 ${\mathbf{D}}_N^{(n)}$ 之外,還必須在線性算符的常數(shù) $c$ 上乘以 $N$ 階單位矩陣 ${\mathbf{I}}_N$, 得到 $N$ 階方陣 $\mathbf{L}$, 即: $\mathbf{L}=a{\mathbf{D}}_N^{(2)}+b{\mathbf{D}}_N+c{\mathbf{I}}_N$ .注意: Matlab 對(duì)矩陣與常數(shù)之間的加法是有特殊定義的, 如果忘記在 $c$ 上乘以 ${\mathbf{I}}_N$ 會(huì)導(dǎo)致錯(cuò)誤, 因?yàn)榫€性算符矩陣 $\mathbf{L}$ 與向量 $\mathbf{u}$ 相乘 $\mathbf{L}\mathbf{u}=\left(a{\mathbf{D}}_N^{(2)}+b{\mathbf{D}}_N+c{\mathbf{I}}_N\right)\mathbf{u}=a{\mathbf{D}}_N^{(2)}\mathbf{u}+b{\mathbf{D}}_N\mathbf{u}+c\mathbf{u}\ne \left(a{\mathbf{D}}_N^{(2)}+b{\mathbf{D}}_N+c\right)\mathbf{u}$ .對(duì)非線性項(xiàng) $N(u)$, 替換 ${?}^n?x^n$ 為 ${\mathbf{D}}_N^{(n)}$ 時(shí), ${\mathbf{D}}_N^{(n)}$ 與向量 $\mathbf{u}$ 之間的乘法是矩陣乘法, 在 Matlab 中用 *表示.而 $u^3$ 應(yīng)理解為向量 $\mathbf{u}$ 中的每個(gè)元素的立方, 這在 Matlab 中用 . ^3 表示.類似地, $f(\mathbf{x})?\left({\mathbf{D}}_N\mathbf{u}\right)$ 也應(yīng)處理為 $f(\mathbf{x})$ 和向量 ${\mathbf{D}}_N\mathbf{u}$ 中的每個(gè)對(duì)應(yīng)元素的相乘, 在 Matlab 中用 .*表示.所以此時(shí) $Lu+N(u)$ 在 Matlab 中應(yīng)寫為:

(a*D2+b*D1+c*I)*u+u.^3+u.^3.*(D2*u)+f.*(D1*u)

換句話說(shuō), 只有譜求導(dǎo)矩陣 ${\mathbf{D}}_N^{(n)}$ 與向量 $\mathbf{u}$ 之間的乘法是矩陣運(yùn)算, 其余的運(yùn)算均是在矩陣或向量的元素之間進(jìn)行的.但也有一個(gè)例外,由于線性算符 $L$ 中的常數(shù) $c$ 在矩陣化時(shí)乘了單位矩陣 ${\mathbf{I}}_N$, 所以實(shí)際上 $c{\mathbf{I}}_N$ 與向量 $\mathbf{u}$ 之間也是矩陣乘法, 這樣才能將其與線性算符 $L$ 中的其他項(xiàng)相加: $\mathbf{L}=a{\mathbf{D}}_N^{(2)}+b{\mathbf{D}}_N+c{\mathbf{I}}_N$ .

若方程式 % \eqref{4-26} 再增加一個(gè)維度, 等號(hào)右邊包含了 $u(x,y,t)$ 的 ${?}^n/?x^n$ 和 ${?}^n/?y^n$, 那么求解步驟中的一些細(xì)節(jié)就需要推廣到二維空間上.在 $x$ 軸、 $y$ 軸上的計(jì)算區(qū)間內(nèi)分別取等間距的 $N$ 個(gè)位置 $\mathbf{x}={\left(x_1,x_2,\dots ,x_N\right)}^{\mathrm{T}}$ 以及 $\mathbf{y}={\left(y_1,y_2,\dots ,y_N\right)}^{\mathrm{T}}$, 于是在 $x?y$ 平面上的計(jì)算區(qū)域內(nèi)就得到了 $N^2$ 個(gè)位置的坐標(biāo) $\left(x_1,y_1\right),\left(x_1,y_2\right),\dots ,\left(x_1,y_N\right),\left(x_2,y_1\right),\dots ,\left(x_2,y_N\right),\dots ,\left(x_N,y_N\right)$ .相應(yīng)地, 函數(shù) $u$ 可以被離散化為 $N$ 階方陣（第一種形式）或 $N^2$ 維列向量（第二種形式）, 如式 % \eqref{4-28} 、式 % \eqref{4-29} 所示:

$$\begin{array}{rl}& {\mathbf{u}}_{N\times N}=\left(\begin{array}{cccc}{u}_{11}& u_{12}& ?& u_{1N}\\ u_{21}& u_{22}& ?& u_{2N}\\ ?& ?& ?& ?\\ u_{N1}& u_{N2}& ?& u_{NN}\end{array}\right)\end{array}$$

$${\mathbf{u}}_{N^2\times \mathbf{l}}={\left(\begin{array}{cccccccccc}{u}_{11},& u_{21},& \dots ,& u_{N1},& u_{12},& u_{22},& \dots ,& u_{N2},& \dots ,& u_{NN}\end{array}\right)}^{\mathrm{T}}$$

式 % \eqref{4-28} 中, 第 $i$ 行 $j$ 列的元素 $u_{ij}$ 對(duì)應(yīng)的坐標(biāo)是 $\left(x_j,y_i\right)$ .列向量 % \eqref{4-29} 的第 1 到第 $N$個(gè)元素對(duì)應(yīng)于式 % \eqref{4-28} 的第 1 列, 第 $N+1$ 到第 $2N$ 個(gè)元素對(duì)應(yīng)于式 % \eqref{4-28} 的第 2 列……依此類推. 另外, 可以通過(guò) Matlab 中的 reshape 函數(shù)在二者間進(jìn)行轉(zhuǎn)換.

收斂階說(shuō)明

收斂階定義

$$order=\{\begin{array}{c}{\log }_2[E(\tau )/E(\tau /2)],\text{?in?time?}\\ {\log }_2[E(N)/E(2N)],\text{?in?space.?}\end{array}$$

時(shí)間方向	空間方向
斜率含義	斜率含義
$$\begin{array}{rl}k& =\frac{{\log }_{2}E(\tau )?{\log }_{2}E\left(\frac{\tau }{2}\right)}{{\log }_2\tau ?{\log }_2\frac{\tau }{2}}\\ & =\frac{{\log }_2\frac{E(\tau )}{E\left(\frac{\tau }{2}\right)}}{{\log }_2\frac{\tau }{\frac{\tau }{2}}}\\ & ={\log }_2\frac{E(\tau )}{E\left(\frac{\tau }{2}\right)}\\ & =\text{?order?}\end{array}$$	$$\begin{array}{rl}k& =\frac{{\log }_{2}E(N)?{\log }_{2}E\left(2N\right)}{N?2N}\\ & =?\frac{1}{N}{\log }_2\frac{E(N)}{E\left(2N\right)}\\ & =?\frac{1}{N}\text{?order?}\end{array}$$
時(shí)間步長(zhǎng)$\tau$與誤差$E(\tau )$的關(guān)系	空間剖分$N$與誤差$E(N)$的關(guān)系
$$\begin{array}{rl}{\log }_{2}E(\tau )& =k{\log }_2\tau +b\\ {\log }_{2}E(\tau )& ={\log }_2{\tau }^k+{\log }_2{2}^b\\ {\log }_{2}E(\tau )& ={\log }_2{2}^b{\tau }^k\\ E(\tau )& =2^b{\tau }^k\end{array}$$	$$\begin{array}{rl}{\log }_{2}E(N)& =kN+b\\ {\log }_{2}E(N)& ={\log }_2{2}^{kN+b}\\ E(N)& =2^{kN+b}\\ E(N)& =2^b{2}^{kN}\end{array}$$