牛頓法──非線性方程的求根方法

本文的閱讀等級:初級

牛頓法 (Newton’s method) 或稱牛頓─拉弗森法 (Newton-Raphson method) 是一個極有效的非線性方程 f(x)=0 的求根方法。令 f:[a,b]\to\mathbb{R} 為一連續可導函數。設 x_0f(x)=0 的一根的估計值,寫出泰勒級數

\displaystyle f(x)=0=f(x_0)+f'(x_0)(x-x_0)+O((x-x_0)^2)

如果 \vert x-x_0\vert 足夠小,我們可以忽略截斷 (truncated) 誤差 O((x-x_0)^2)。解線性方程 f(x_0)+f'(x_0)(x-x_0)=0 可得近似根:

\displaystyle x_1=x_0-\frac{f(x_0)}{f'(x_0)}

上式的幾何意義是 (x_1,0) 為函數 f 於點 (x_0,f(x_0)) 的切線與x-軸的交點。我們期待 x_1x_0 更接近真實根,遂以迭代程序連續逼近:對於 k=0,1,2,\ldots

\displaystyle x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}

 
下圖顯示牛頓法的迭代過程:藍線是給定的函數 f,紅線是切線。設定目標誤差 \epsilon>0,持續迭代直到 \vert x_{k+1}-x_k\vert<\epsilon。明顯地,若 f'(x_k)=0,牛頓法將被迫中止。往下閱讀前,讀者不妨利用線上牛頓法計算器來解方程式,並嘗試更換不同的初始值以體會牛頓法的收斂性能。多數的時候,牛頓法都能成功地收斂至一根。不過也有少數例外,譬如,

\displaystyle x^3-x+\frac{\sqrt{2}}{2}=0

迭代公式為

\displaystyle x_{k+1}=x_k-\frac{x_k^3-x_k+\sqrt{2}/2}{3x_k^2-1}

x_0=0,則 x_1=\sqrt{2}/2x_2=0\{x_k\} 數列將來回震盪,無法收斂。若 x_0=-1,經過7次迭代便收斂至 -1.25107862158365

Newton Iteration
Image source: Wikipedia, by Ralf Pfeifer

 
牛頓法是否能收斂取決於初始值與根的距離遠近,下面的收斂定理說明初始值的重要性。

收斂性定理:令 f 是定義於 [a,b] 的二次可導函數。若 p\in[a,b] 使得 f(p)=0f'(p)\neq 0,則存在一區間 I=[p-\delta,p+\delta]\delta>0,使得對於任一 x_0\in I,牛頓法產生的數列 \{x_k\} 收斂至 p

我們要證明當 k\to\inftyx_k\to p。令

\displaystyle g(x)=x-\frac{f(x)}{f'(x)}

牛頓法的算式 x_{k+1}=g(x_k) 其實是一種不動點迭代 (fixed point iteration)。因為 f(p)=0,即知 g(p)=p,故 pg 的一個不動點。因為 f'(p)\neq 0f' 是連續函數,存在 \epsilon>0 使每一 x\in [p-\epsilon,p+\epsilon]\subset[a,b] 都有 f'(x)\neq 0。這說明 g(x) 在區間 [p-\epsilon,p+\epsilon] 定義良好且連續。所以,

\displaystyle g'(x)=1-\frac{f'(x)f'(x)-f(x)f''(x)}{(f'(x))^2}=\frac{f(x)f''(x)}{(f'(x))^2}

即有 g'(p)=0。任意挑選正數 \rho<1,根據連續性,必存在區間 I=[p-\delta,p+\delta] 使得 \vert g'(x)\vert\le\rhox\in I。對於 x\in I,中值定理 (mean value theorem) 指出存在 c\in I 滿足

\displaystyle \vert g'(c)\vert=\frac{\vert g(p)-g(x)\vert}{\vert p-x\vert}

\displaystyle \vert p-g(x)\vert=\vert g(p)-g(x)\vert=\vert p-x\vert \vert g'(c)\vert\le\delta\rho<\delta

上式可寫成 p-\delta<g(x)<p+\delta,因此 g(x)\in I。最後套用不動點定理[1],證明 x_{k+1}=g(x_k) 收斂至 p

 
透過誤差分析可以推算出牛頓法的收斂速度。設 pf(x)=0 的一根,f(p)x_k 的泰勒公式為

\displaystyle f(p)=f(x_k)+f'(x_k)(p-x_k)+\frac{1}{2}f''(\xi_k)(p-x_k)^2

其中 \xi_k 位於 x_kp 之間。因為 f(p)=0,上式等號兩邊同除以 f'(x_k),可得

\displaystyle 0=\frac{f(x_k)}{f'(x_k)}+(p-x_k)+\frac{f''(\xi_k)}{2f'(x_k)}(p-x_k)^2

又因為 x_{k+1}=x_k-f(x_k)/f'(x_k)

\displaystyle p-x_{k+1}=-\frac{f''(\xi_k)}{2f'(x_k)}(p-x_k)^2

e_k=p-x_k。代入上式,等號兩邊取絕對值即有

\displaystyle \vert e_{k+1}\vert =\frac{\vert f''(\xi_k)\vert}{2\vert f'(x_k)\vert}e_k^2

故知牛頓法具有平方收斂的性能。粗略地說,每迭代一次,新估計值的有效數字將增加一倍。

 
接著我們將單變量函數推廣至多變量函數。令 F:D\to\mathbb{R}^n 為一可導向量函數,其中 D\subset\mathbb{R}^n。對於 \mathbf{x}=(x_1,\ldots,x_n)^T\in D

\displaystyle F(\mathbf{x})=\begin{bmatrix} f_1(x_1,\ldots,x_n)\\ \vdots\\ f_n(x_1,\ldots,x_n) \end{bmatrix}

\mathbf{x}_k=(x_{k1},\ldots,x_{kn})^TF(\mathbf{x})=\mathbf{0} 的一根的估計值。考慮 f_i(\mathbf{x})\mathbf{x}_k 的泰勒級數:

\displaystyle f_i(\mathbf{x})=0=f_i(\mathbf{x}_k)+\sum_{j=1}^n\frac{\partial f_i}{\partial x_j}(\mathbf{x}_k)(x_j-x_{kj})+\cdots,~~~i=1,\ldots,n

如果忽略二次以上所有項,上面 n 個線性方程可用矩陣表示如下:

\displaystyle J(\mathbf{x}_k)(\mathbf{x}-\mathbf{x}_k)=-F(\mathbf{x}_k)

其中 J(\mathbf{x}_k) 是 Jocobian 矩陣 (見“Jacobian 矩陣與行列式”)

J(\mathbf{x}_k)=\begin{bmatrix} \frac{\partial f_1}{\partial x_1}(\mathbf{x}_k)&\frac{\partial f_1}{\partial x_2}(\mathbf{x}_k)&\cdots&\frac{\partial f_1}{\partial x_n}(\mathbf{x}_k)\\[0.5em] \frac{\partial f_2}{\partial x_1}(\mathbf{x}_k)&\frac{\partial f_2}{\partial x_2}(\mathbf{x}_k)&\cdots&\frac{\partial f_2}{\partial x_n}(\mathbf{x}_k)\\ \vdots&\vdots&\ddots&\vdots\\ \frac{\partial f_n}{\partial x_1}(\mathbf{x}_k)&\frac{\partial f_n}{\partial x_2}(\mathbf{x}_k)&\cdots&\frac{\partial f_n}{\partial x_n}(\mathbf{x}_k) \end{bmatrix}

J(\mathbf{x}_k) 是可逆的,牛頓法的迭代公式為

\displaystyle \mathbf{x}_{k+1}=\mathbf{x}_k-J(\mathbf{x}_k)^{-1}F(\mathbf{x}_k)

在實際應用時,你不需要真的算出逆矩陣 J(\mathbf{x}_k)^{-1},使用高斯消去法解出線性方程 J(\mathbf{x}_{k})\mathbf{z}_k=-F(\mathbf{x}_k),再代入 \mathbf{x}_{k+1}=\mathbf{x}_k+\mathbf{z}_k 即可 (見“高斯消去法”)。

 
舉一個例子[2]

F(x_1,x_2)=\begin{bmatrix} x_1\\ \frac{10x_1}{x_1+0.1}+2x_2^2 \end{bmatrix}=\mathbf{0}

有唯一根 (x_1,x_2)=(0,0)^T。Jacobian 矩陣為

J(x_1,x_2)=\begin{bmatrix} 1&0\\ \frac{1}{(x_1+0.1)^2}&4x_2 \end{bmatrix}

x_2=0,則 J(x_1,x_2) 是不可逆的。若 \mathbf{x}_k=(0,y)^T,則

\mathbf{x}_{k+1}=\begin{bmatrix} 0\\ y \end{bmatrix}-\begin{bmatrix} 1&0\\ 100&4y \end{bmatrix}^{-1}\begin{bmatrix} 0\\ 2y^2 \end{bmatrix}=\begin{bmatrix} 0\\ \frac{y}{2} \end{bmatrix}

因此,\Vert\mathbf{x}_{k+1}\Vert=\frac{1}{2}\Vert\mathbf{x}_k\Vert,此例牛頓法具有線性的收斂性。

 
註解
[1] 不動點定理 (fixed point theorem):設 g 是定義於 [a,b] 的一次可導函數。若對於任一 x\in[a,b]g(x)\in[a,b]\vert g'(x)\vert<1,則 x_{k+1}=g(x_k) 收斂至 p,其中 p\in[a,b]g 的一個不動點。

[2] 取自 M.J.D. Powell, A Hybrid Method for Non-Linear Equations, in P. Rabinowitz (ed): Numerical Methods for Non-Linear Algebraic Equations, 1970, pp 87.

廣告
本篇發表於 特別主題 並標籤為 , , , 。將永久鏈結加入書籤。

8 Responses to 牛頓法──非線性方程的求根方法

  1. yufenglyu 說道:

    周老师,正文最后一行“使用高斯消去法解出线性方程 J(\mathbf{x}_{k+1})\mathbf{z}_k=-F(\mathbf{x}_k)”,是不是应该是 J(\mathbf{x}_{k})\mathbf{z}_k=-F(\mathbf{x}_k)

  2. 董瀚林 說道:

    周老师您好,很荣幸拜读您的文章。但有一处希望与您探讨:您在本文收敛性证明处,对于迭代公式的表述写到g(x)=x-f(x)/f'(x),此处g为f的反函数,那么表述成g(y)是否更加妥当?不对之处恳请您批评指正。

  3. 查理 說道:

    請教一下:為何牛頓法n個聯立方程式時,會碰到無解的請況呢?是跟Jacobian有關?

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s