線性代數在數值分析的應用 (一):兩點邊值問題

本文的閱讀等級:初級

數值分析 (numerical analysis) 是一門研究連續數學問題演算法的科目,內容包羅萬象,如解非線性方程、解線性或非線性聯立方程組、插值 (interpolation) 與曲線配適 (curve fitting)、數值微分與積分、常微分方程與偏微分方程的數值解、邊值 (boundary-value) 與特徵值 (characteristic-value) 問題等。線性代數和數值分析有極密切的關係,一方面數值分析涵蓋一些數值線性代數問題,另一方面數值分析的部分演算法建立於線性代數之上。本文介紹線性代數在數值分析中的一個簡單應用──將兩點邊值 (two-point boundary-value) 問題轉換為線性聯立方程組的求解問題。

 
考慮下面的二階微分方程:

{x}''-x=t

兩邊值為 x(0)=2x(1)=-1。我們的目標是計算 x 於區間 [0,1]n 個等距間隔點 t_i 的近似值,t_i=i/(n+1)i=1,2,\ldots,n

 
求兩點邊值問題數值解的第一步驟是將連續函數的微分以差分取代。設函數所定義的區間為 [a,b],將此區間分割為 n+1 個等份,令 h=\frac{b-a}{n+1},切割的等間隔點為 t_i=a+ihi=0,1,\ldots,n+1。考慮在 t_i+ht_i-h 的泰勒展開式 (Taylor series expansion):

\begin{aligned}\displaystyle  x(t_i+h)&=x(t_i)+x'(t_i)h+ \frac{h^2}{2!}x'' (t_i)+\frac{h^3}{3!}x''' (t_i)+\cdots\\    x(t_i-h)&=x(t_i)-x'(t_i)h+ \frac{h^2}{2!}x'' (t_i)-\frac{h^3}{3!}x''' (t_i)+\cdots\end{aligned}

將上面二式相減與相加,分別得到

\begin{aligned}\displaystyle  x'(t_i)&=\frac{x(t_i+h)-x(t_i-h)}{2h}+O(h^2)\\    x''(t_i)&=\frac{x(t_i+h)-2x(t_i)+x(t_i-h)}{h^2}+O(h^2)\end{aligned}

其中符號 O(h^k) 代表包含 k 次冪和更高次冪的項。為簡化書寫,令 x_{i}=x(t_i)。如果捨去 O(h^2),結果稱為中央差分近似 (central difference approximation),如下:

\begin{aligned}\displaystyle  x'(t_i)&\approx\frac{x_{i+1}-x_{i-1}}{2h}\\    x''(t_i)&\approx\frac{x_{i+1}-2x_i+x_{i-1}}{h^2}\end{aligned}

已知 a=0b=1,故 t_i=ihh=1/(n+1)。當 t=t_i 時,將 x''(t_i) 的近似式代入原方程式:

\displaystyle  \frac{x_{i+1}-2x_i+x_{i-1}}{h^2}-x_i=t_i

整理為

x_{i-1}-(2+h^2)x_i+x_{i+1}=h^2t_i=ih^3,~~i=1,2,\ldots,n

另外還需考慮已知邊值 x_0=x(t_0)=x(0)=2x_{n+1}=x(t_{n+1})=x(1)=-1

 
設未知向量 \mathbf{x}=(x_1,x_2,\ldots,x_n)^T,上述 n 個聯立方程式可表示為矩陣形式 A\mathbf{x}=\mathbf{b}An\times n 階,如下:

\begin{aligned}  &\begin{bmatrix}    -2-h^2&1&0&\cdots&0&0&0\\    1&-2-h^2&1&\cdots&0&0&0\\    0&1&-2-h^2&\cdots&0&0&0\\    \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\    0&0&0&\cdots&-2-h^2&1&0\\    0&0&0&\cdots&1&-2-h^2&1\\    0&0&0&\cdots&0&1&-2-h^2    \end{bmatrix}\begin{bmatrix}  x_1\\  x_2\\  x_3\\  \vdots\\  x_{n-2}\\  x_{n-1}\\  x_n  \end{bmatrix}\\  &=\begin{bmatrix}    h^3-2\\    2h^3\\    3h^3\\    \vdots\\    (n-2)h^3\\    (n-1)h^3\\    nh^3+1    \end{bmatrix}\end{aligned}

解出此線性方程即可得到 x(t) 的近似值,如果增加區間切割數 n+1h 值將隨之減小,可得到更緊密且精確的結果。當然,代價是必須解出一個尺寸較大的 n 階線性方程。觀察上面的係數矩陣形式發現它是一個三對角矩陣,以高斯消去法化簡三對角矩陣僅需消耗 O(n) 計算量 (見“特殊矩陣 (11):三對角矩陣”)。

廣告
本篇發表於 線性代數專欄, 應用之道 並標籤為 , , , , 。將永久鏈結加入書籤。

發表迴響

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

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s