線性代數在數值分析的應用 (二):Dirichlet 問題

本文的閱讀等級:中級

在物理學中,等方向均勻介質的一點的溫度變化由熱傳導方程 (heat equation) 所描述[1]

\displaystyle  \frac{\partial u}{\partial t}=\alpha\nabla^2u

其中 u=u(x,y,z,t) 是點 (x,y,z) 於時間 t 的溫度,\alpha 是一正數,稱為熱擴散率 (thermal diffusivity),\nabla^2 是 Laplace 算子 (或稱 Laplacian),定義如下:

\displaystyle  \nabla^2u=\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}

淺白地說,Laplace 算子度量一點的函數值與其鄰近點的差異。若點 (x,y,z) 處於穩態,即該點溫度不隨時間改變,則 u=u(x,y,z) 滿足 Laplace 方程 \nabla^2u=0。如果二階連續可導函數 u:S\to\mathbb{R} (S\subset\mathbb{R}^3 為一開集) 滿足 Laplace 方程,我們稱之為調和函數 (harmonic function)。本文將探討簡化後的二維 Dirichlet 問題[2]:尋找一調和函數 u(x,y),使其在一正方形區域內所有點皆滿足 \nabla^2u=\partial^2u/\partial x^2+\partial^2u/\partial y^2=0,且邊界滿足給定條件 u(x,y)=c(x,y)

 
我們採用數值分析方法,首先將 Laplace 算子予以離散化。令正方形區域等分切割為 (n+1)\times(n+1) 個小正方形,h 表示小正方形的邊長。設離散點的座標為 (x_i,y_j)i,j=0,1,\ldots,n+1,其中內點的指標為 i,j=1,2,\ldots,n,邊界點的 i,j 指標至少有一個等於 0n+1。下圖顯示 n=3i,j 指標分布:

Dirichlet 問題1

 
在內點 (x_i,y_j),我們以中央差分近似 (central difference approximation) 取代 \partial^2u/\partial x^2\partial^2u/\partial y^2 (見“線性代數在數值分析的應用 (一):兩點邊值問題”):

\displaystyle\begin{aligned}  \left.\frac{\partial^2u}{\partial x^2}\right|_{(x_i,y_j)}&=\frac{u(x_i-h,y_j)-2u(x_i,y_j)+u(x_i+h,y_j)}{h^2}+O(h^2)\\  \left.\frac{\partial^2u}{\partial y^2}\right|_{(x_i,y_j)}&=\frac{u(x_i,y_j-h)-2u(x_i,y_j)+u(x_i,y_j+h)}{h^2}+O(h^2).\end{aligned}

將上面兩式相加,通乘 h^2,並令內點 (x_i,y_j) 滿足調和函數條件 \partial^2u/\partial x^2+\partial^2u/\partial y^2=0,可得

\displaystyle  4u_{ij}=u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}+O(h^4),~~~i,j=1,2,\ldots,n

上式中,我們採用了簡化符號 u_{ij}=u(x_i,y_j)。如果捨去 O(h^4),則得下列5-點差分 Laplace 方程組:

\displaystyle  4u_{ij}-(u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1})=0,~~~i,j=1,2,\ldots,n

下圖顯示5個點的相對位置:

Dirichlet 問題2

離散 Laplace 算子可以視為以內點 (x_i,y_j) 為中心與下列 3\times 3 階矩陣的摺積 (convolution):

\displaystyle  \left[\!\!\begin{array}{rrr}  0&-1&0\\  -1&4&-1\\  0&-1&0  \end{array}\!\!\right]

(i,j) 為邊界指標,則 u_{i,j}=c_{i,j},其中 c_{i,j}=c(x_i,y_j)。我們可以將 Laplace 方程組以矩陣表示為 L\mathbf{u}=\mathbf{b},其中 L 是一 n^2\times n^2 矩陣,\mathbf{u}\mathbf{b}n^2 維向量。以 n=3 為例,

\displaystyle  \left[\!\!\begin{array}{rrrcrrrcrrr}  4&-1&0&\vline&-1&0&0&\vline&0&0&0\\  -1&4&-1&\vline&0&-1&0&\vline&0&0&0\\  0&-1&4&\vline&0&0&-1&\vline&0&0&0\\ \hline  -1&0&0&\vline&4&-1&0&\vline&-1&0&0\\  0&-1&0&\vline&-1&4&-1&\vline&0&-1&0\\  0&0&-1&\vline&0&-1&4&\vline&0&0&-1\\\hline  0&0&0&\vline&-1&0&0&\vline&4&-1&0\\  0&0&0&\vline&0&-1&0&\vline&-1&4&-1\\  0&0&0&\vline&0&0&-1&\vline&0&-1&4  \end{array}\!\!\right]\begin{bmatrix}  u_{11}\\  u_{12}\\  u_{13}\\ \hline  u_{21}\\  u_{22}\\  u_{23}\\ \hline  u_{31}\\  u_{32}\\  u_{33}  \end{bmatrix}=\begin{bmatrix}  c_{01}+c_{10}\\  c_{02}\\  c_{03}+c_{14}\\ \hline  c_{20}\\  0\\  c_{24}\\ \hline  c_{30}+c_{41}\\  c_{42}\\  c_{43}+c_{34}  \end{bmatrix}

 
係數矩陣 L 稱為離散 Laplacian,具有下列對稱三對角分塊形式:

\displaystyle  L=\begin{bmatrix}  A&-I&&&\\  -I&A&-I&&\\  &\ddots&\ddots&\ddots&\\  &&-I&A&-I\\  &&&-I&A  \end{bmatrix}

其中 An\times n 階對稱三對角矩陣

\displaystyle  A=\left[\!\!\begin{array}{rrcrr}  4&-1&&&\\  -1&4&-1&&\\  &\ddots&\ddots&\ddots&\\  &&-1&4&-1\\  &&&-1&4  \end{array}\!\!\right]

離散 Laplacian L 是一正定矩陣。下面我們直接計算證明 L 的特徵值皆為正數 (見“正定矩陣的性質與判別方法”)。根據三對角矩陣的特徵值公式 (見“每週問題 February 15, 2010”),A 的特徵值為

\displaystyle  \lambda_i=4-2\cos\left(\frac{i\pi}{n+1}\right),~~~i=1,2,\ldots,n

另一方面,A 是一實對稱矩陣,故可正交對角化為 Q^TAQ=D=\hbox{diag}(\lambda_1,\ldots,\lambda_n),其中 QQ^T=Q^TQ=I。令 Un^2\times n^2 階分塊主對角正交矩陣

\displaystyle  U=\begin{bmatrix}  Q&0&\cdots&0\\  0&Q&\cdots&0\\  \vdots&\vdots&\ddots&\vdots\\  0&0&\cdots&Q  \end{bmatrix}

不難驗證下列相似變換:

\displaystyle  U^TLU=\begin{bmatrix}  D&-I&&&\\  -I&D&-I&&\\  &\ddots&\ddots&\ddots&\\  &&-I&D&-I\\  &&&-I&D  \end{bmatrix}=B

接著利用排列矩陣置換行列指標。設 P 為一 n^2\times n^2 階排列矩陣使得 B 相似於一分塊對角矩陣,如下:

\displaystyle  P^TBP=\begin{bmatrix}  M_1&0&\cdots&0\\  0&M_2&\cdots&0\\  \vdots&\vdots&\ddots&\vdots\\  0&0&\cdots&M_n  \end{bmatrix}=M

其中 n\times n 階主對角分塊 M_i 為對稱三對角矩陣

\displaystyle  M_i=\begin{bmatrix}  \lambda_i&-1&&&\\  -1&\lambda_i&-1&&\\  &\ddots&\ddots&\ddots&\\  &&-1&\lambda_i&-1\\  &&&-1&\lambda_i  \end{bmatrix},~~~i=1,2,\ldots,n

n=3 為例,P 對應排列

\pi=\begin{pmatrix}  1&2&3&4&5&6&7&8&9\\  1&4&7&2&5&8&3&6&9  \end{pmatrix}

亦即 (見“特殊矩陣 (16):排列矩陣”)

\displaystyle  P=\begin{bmatrix}  1&0&0&0&0&0&0&0&0\\  0&0&0&1&0&0&0&0&0\\  0&0&0&0&0&0&1&0&0\\  0&1&0&0&0&0&0&0&0\\  0&0&0&0&1&0&0&0&0\\  0&0&0&0&0&0&0&1&0\\  0&0&1&0&0&0&0&0&0\\  0&0&0&0&0&1&0&0&0\\  0&0&0&0&0&0&0&0&1  \end{bmatrix}

因為 L 相似於 BB 相似於 M,故 LM 有相同的特徵值。分塊對角矩陣 M 的特徵值由所有分塊 M_i 的特徵值構成。再次使用三對角矩陣的特徵值公式,M_i 的特徵值為 \lambda_i-2\cos(j\pi/(n+1))j=1,\ldots,n。代入 \lambda_i 的表達式,即得 M (也是 L) 的特徵值

\displaystyle  \lambda_{ij}=4-2\left(\cos\left(\frac{i\pi}{n+1}\right)+\cos\left(\frac{j\pi}{n+1}\right)\right),~~~i,j=1,2,\ldots,n

使用半角公式 1-\cos\theta=2\sin^2(\theta/2),離散 Laplacian L 的特徵值如下:

\displaystyle  \lambda_{ij}=4\left(\sin^2\left(\frac{i\pi}{2n+2}\right)+\sin^2\left(\frac{j\pi}{2n+2}\right)\right),~~~i,j=1,2,\ldots,n

因為每一 \lambda_{ij}>0,證明 L 是一正定矩陣。這個結果同時表明 L 可逆,故 \mathbf{u}=L^{-1}\mathbf{b},也就是說,不論邊界條件為何,Dirichlet 問題總是存在唯一解。

 
參考來源:
[1] 維基百科:Heat eqaution
[2] Carl D. Meyer 所著 Matrix Analysis and Applied Linear Algebra,2001,頁 563-565。

This entry was posted in 線性代數專欄, 應用之道 and tagged , , , , , , , . Bookmark the permalink.

2 則回應給 線性代數在數值分析的應用 (二):Dirichlet 問題

  1. wonderlandtommy 說:

    版主好,
    我正在操作三維的離散Laplacian,請問求三維Laplacian的特徵值時應該將三階張量降維矩陣嗎?謝謝版主

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 變更 )

Twitter picture

You are commenting using your Twitter account. Log Out / 變更 )

Facebook照片

You are commenting using your Facebook account. Log Out / 變更 )

Google+ photo

You are commenting using your Google+ account. Log Out / 變更 )

連結到 %s