【機器學習】數值分析02——任意方程求根

任意方程求根

全文目錄

（部落格園）機器學習

（Github）MachineLearning Math

1.簡介

方程和函式是代數數學中最為重要的內容之一，從初中直到大學，我們都在研究著方程與函式，甚至我們將圖形代數化，從而發展出了代數幾何、解析幾何的內容。而在方程與函式中，我們研究其性質最多的，往往就是方程的根（零點），即使是研究方程的極值點、鞍點等，我們無非也只是研究其微商的零點。
我們在初等數學中已經學過許多簡單初等函式、線性方程的求解方法，在本文中，我們重點討論任意方程，尤其是計算困難的非線性方程的求根方法。

2.方程

2.1分類和介紹

方程就是指含有未知數的等式。是表示兩個數學式（如兩個數、函式、量、運算）之間相等關係的一種等式，使等式成立的未知數的值稱為“解”或“根”。在這裡，根據一些性質的不同，我們將方程分成以下幾類：

• 單個方程
  • 線性方程：本質是等式兩邊乘以任何相同的非零數，方程的本質都不受影響。通常認為只含有一次項的方程。
  • 非線性方程：是因變數與自變數之間的關係不是線性的關係的方程。
    • 多項式方程
    • 超越方程：指含有未知量的超越式（指數、對數、三角函式、反三角函式等）的方程。換言之，超越方程中都有無法用含有未知數的多項式、分式或開方表示的式子。
• 多個方程
  • 線性方程組
  • 非線性方程組

2.2方程的零點（根、解）

若有一個值或一些值能夠使得方程 \(f(x)=0\) 成立，那麼這個值就被成為方程的解，也常常被叫做零點和根。

若方程有且只有一個解\(x^*\)，那麼我們稱方程有單根\(x^*\)。

若對於方程\(f(x)=0\)，有\(f(x^*) = 0,f^{'}(x^*)=f^{''}(x^*)=\cdots=f^{(k)}(x^*)=0,f^{(k+1)}(x^*)\neq0\),那麼稱\(x^*\)為方程的k+1重根

PS：若方程是簡單冪函式多項式組成，那麼方程的解的數量應和最高此項的數值一致，因為存在虛根。

3.求根方法

求根的方法基本上大同小異，都是通過區間去逼近方程的根的點。

首先我們說一個定理1：對於實函式方程\(f(x)=0\)，當\(x\in(a,b)\)，且\(f(x)\)在\(x\in(a,b)\)時單調且連續，若\(f(a)\cdot f(b)<0\)，則方程在\(x\in(a,b)\)有且只有一個根。

3.1二分法

3.1.1普通二分法的原理及操作

二分法和演算法中的二分搜尋法非常的類似，取定有根區間\([a,b]\)進行對分，求得\(mid = \frac{a+b}{2}\)進行判斷含根區間，如果\(f(a)\cdot f(mid)<0\),則令\(b=mid\)；反之若\(f(b)\cdot f(mid)<0\),則令\(a=mid\)。當\(|b_n-a_n|<\epsilon\)停止計算返回結果。

產生的截斷誤差為\(|e_{n-1}| = x_{n+1} - x^*\leq[b_n - a_n] = \frac{b_0 - a_0}{2^n}\)。

可以計算出最小迭代次數為\(n = \frac{lg(c_0-a_0)-lg\epsilon}{lg2}\)

程式碼實現（更多語言的程式碼見倉庫中Code資料夾）：

private static double epsilon = 0.001;
// func為函式，寫法如x=>x*x+2*x-1，a，b必須為有效的含根開區間
public static double Binary(Func<double, double> func, double a, double b)
{
    var f1 = func.Invoke(a);
    var f2 = func.Invoke(b);
    if (f1 * f2 > 0)
        throw new ArgumentException("此區間無根或根不唯一");
    double mid = (a + b) / (double)2;
    var fm = func.Invoke(mid);
    if (fm == 0)
        return fm;
    if (f1 * fm < 0)
        b = mid;
    else if (f2 * fm < 0)
        a = mid;
    if (Math.Abs(b - a) <= epsilon)
        return (a + b) / (double)2;
    return Binary(func, a, b);
}

3.1.2普通二分法準確度及速度分析

假設\(\left[a,b\right]\)是二分法的初始區間，在進行n此二分之後，得到的最終根的分佈區間\([a_n,b_n]\)的長度為\((b-a)/2^n\)，我們將其中點作為根的最優估計值，與真實值之間的誤差不超過區間長度的一半。

\[e = \left|x_c-r\right|<\frac{b-a}{2^{n+1}} \]

若\(e<\frac{1}{2}\times10^{-p}\),則精確到小數點後p位，

事實上我們也可以簡單的計算出函式執行的次數為n+2次。

3.2淺談迭代

3.2.1 迭代是什麼

迭代是重複反饋過程的活動，其目的通常是為了逼近所需目標或結果。每一次對過程的重複稱為一次“迭代”，而每一次迭代得到的結果會作為下一次迭代的初始值。
重複執行一系列運算步驟，從前面的量依次求出後面的量的過程。此過程的每一次結果，都是由對前一次所得結果施行相同的運算步驟得到的。這篇文章中對於迭代有一個非常不錯的概述究竟什麼是迭代？

3.2.2 不動點的定義

各位可以嘗試以下操作，

1. 隨意輸入一個數字\(\lambda\)
2. 然後對其進行\(cos\lambda\)運算，
3. 將運算結果作為新的值傳回\(cosx\)函式之中

當你一直重複以上三個操作，你會發現數字最後會定格在0.73908513左右不動了。這是一個非常有趣的現象，事實上我們的這個操作就是\(x_n = cosx_{n-1}\)，而最後不變化的數值實際上就是\(x = cosx\)這個方程的解。

這就是我們不動點的一種實際情況，不動點原理是數學上一個重要的原理，也叫壓縮映像原理或巴拿赫（Banach）不動點定理，完整的表達：完備的度量空間上，到自身的一個壓縮對映存在唯一的不動點。用初等數學可以這麼理解：連續對映\(f\)的定義域包含值域，則存在一個\(x\)使得\(f(x)=x\)。

若某函式滿足\(f(\lambda)=\lambda,\lambda \in R\)，我們就稱\(\lambda\)為函式的一階不動點。同樣的，我們推廣一下，若\(f(f(\lambda))=\lambda,\lambda \in R\)，則稱為二階不動點，一階不動點必定是二階不動點。

不動點的存在定理：若某函式\(y=f(x)\)與\(y=x\)存在至少一個交點，那麼函式必然存在不動點。

3.2.3 函式的相似性

我們常常說圖形之間的相似，事實上函式也有其相似的定義。若函式\(f(x)\)對其換元，令\(t = \varphi(x),x=\varphi^{-1}(t),f(x)=>f(\varphi^{-1}(t))\)。我們的不動點是一個基於迭代的過程，迭代就是讓\(f(x_0) = x_1 = t,f(x_1)=f(f(x_0))=f(t)=x_2...\)的過程，如果我們換個角度去思考，迭代不就是和我們剛剛提到的換元是一個道理嗎？那麼我們思考一下，如果對於二次迭代\(f(f(x))\)，如果說我們把函式完全看成是一個複合函式，那麼我們現在需要做的就是試圖將函式寫成\(f(\varphi^{-1}(t))\)的形式。

現在做出如下變換，將外層的\(f(x)\)換元為\(f(\varphi^{-1}(t))\)，那麼二次迭代式就變成了\(g(t) = \varphi(f(\varphi^{-1}(t)))\)。隨後對我們的\(g(x)\)再進行迭代就會得到\(g(g(t)) = \varphi(f(f(\varphi^{-1}(t)))\)，利用數學歸納法計算n次迭代式，那麼就會得到\(g_n(t) = \varphi(f_n(\varphi^{-1}(t)))\)。

我們在這給出總結和定義：若函式\(g(x),f(x),\varphi(x)\)，若有\(g(x)=\varphi(f(\varphi^{-1}(x)))\)，那麼稱函式\(g(x)\)和\(f(x)\)通過\(\varphi(x)\)相似，記作\(f\sim g\)，\(\varphi(x)\)稱為相似的橋函式，且之前我們的過程得出了一個重要的結論就是\(g\sim f=>g^n\sim f^n\)，同時，相似函式的不動點是完全一致的，證明如下：

\[若函式g(x)的不動點為x_0，則有 \]

\[g(x_0)=x_0=\varphi(f(\varphi^{-1}(x_0))) \]

\[根據反函式性質有\varphi(\varphi^{-1}(x))=x \]

\[故f(\varphi^{-1}(x_0)) = \varphi^{-1}(x_0)=x_0,\varphi^{-1}(x_0)是f(x)的一個不動點 \]

通過上述得出的不動點，很容易通過剛才講的內容得出\(\varphi^{-1}(x_0)\)實際上也是\(g(x)\)換元后得到的不動點的位置，因此相似函式的不動點是一一對應的。

3.2.4 迭代收斂速度

迭代是一種逼近、利用數列、函式收斂的性質進行求解的方法，那麼對於收斂的速度，我們很難直接從數值看出速度，因為收斂過程中，誤差的變化是越來越小的，因此我們再次進行一個比階。

\[\exists p\in N \And C>0,\displaystyle \lim_{k \to \infty}\frac{e_{k+1}}{e_k^p} = C \]

其中：

\[e_k = |x_k-x^*| \]

\[p= \begin{cases} p=1,線性收斂\\ p=2,平方收斂\\ p>2,高階收斂 \end{cases} \]

並且有一個定理：若\(\varphi(x)\)在所求的根\(x^*\)鄰域有連續的二階導數，並且有\(|\varphi^{'}(x)|<1\)，有以下特點：

1. 當\(\varphi^{'}(x)\neq0\)，那麼迭代過程線性收斂

2. 當\(\varphi^{'}(x)=0,\varphi^{''}(x)\neq0\)，那麼迭代是平方收斂的。

證明過程留給讀者，只需要利用泰勒展開是便可以證明該定理。

3.3不動點迭代法

3.3.1不動點迭代法的原理及操作

迭代法是將方程求根問題轉換成了函式求交點問題再轉換成數列求極限的問題，這個過程中，對方程進行一個巧妙的變換之後，方程就可以在若干次迭代之後解出一個近似解。

操作方法如下：首先將方程\(f(x)=0\)改寫成\(x = g(x)\)的形式，這樣就可以將方程的解看成是函式\(y=x\)和\(y=g(x)\)的交點。給定初始值\(x_0\)後，則\(x_1 = g(x_0)\)，不斷重複這個過程，若\(\displaystyle \lim_{k \to \infty}x_k\)存在，則迭代便可以達到使得\(x_k\)趨於交點，而這個交點就是我們剛才講了那麼久的不動點。

我們即使構造出了\(x=g(x)\)的迭代式，往往由於我們的變換隻是一些簡單的移項、通分等四則運算獲得，使得迭代式也並不是一個很容易被計算的函式，這個時候，我們之前講到的那麼多函式相似就派上了用場，我們可以利用函式相似性去尋找一個更為簡單的函式\(\varphi(x)\sim g(x)\)，由於相似的性質，不動點並不會發生變化。不過，我們需要保證最終的迭代函式在指定的含根區間記憶體在一階導數且導數值\(|\varphi^{'}(x)|<1\)，否則迭代函式將會不收斂。

迭代過程收斂的情況如下圖所示：

如果說我們的導數不滿足條件，那麼我們的迭代將會發散：

程式碼樣例

// func是迭代函式而不是實際函式
public static double Iterative(Func<double, double> func, double x)
{
    double temp = func.Invoke(x);
    if (Math.Abs(temp - x) <= epsilon)
    {
        return temp;
    }
    x = temp;
    return Iterative(func, x);
}

3.3.2迭代法的收斂性證明

在這裡，我們將證明迭代法求根的合理性和可行性。

前提條件：設函式\(\varphi(x), x\in[a,b]\)有連續的一階導數，並且滿足以下條件：

• \(\forall x\in [a,b],\varphi(x)\in[a,b]\)

• \(\exists L \in (0,1),\forall x\in[a,b],|\varphi^{'}(x)|\leq L\)

要證明和解決以下命題和問題：

• 在\(x\in[a,b],\exists x^*,\varphi(x^*) = 0\)

• \(\forall x_0\in[a,b]\)，迭代過程\(x_{k+1} = \varphi(x_k)\)均收斂與\(x^*\)

• 求解誤差分析式

現在開始證明

1：證明在區間內有且只有一個根存在：

證：在\(x \in [a,b]\)時，\(\varphi^{'}(x)\)存在，所以有\(\varphi(x)\)連續，於是可以作\(g(x) = x - \varphi(x)\)，易知\(g(x)\)連續。

又因為\(\varphi(x) \in [a,b]\)，且\(g(a)*g(b)<0\),故存在實根，使得\(x=\varphi(x)\)。

利用反證法：若在[a,b]上還有一實根\(\bar{x}\)，那麼通過拉格朗日中值定理必定有：

\[x^*-\bar{x} = \varphi(x^*)-\varphi(\bar{x}) = \varphi^{'}(\xi)(x^*-\bar{x})\Longrightarrow \varphi^{'}(\xi) = 1 \]

顯然與假設不符合。

2：證明這個根收斂於\(x^*\)

根據拉格朗日中值定理，容易知道

\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]

又因\(|\varphi^{'}(x)|\leq L\)，故易得：

\[|x^*-x_{k+1}|\leq|L(x^*-x_k)|=|L^2(x^*-x_{k-1})=\cdots=|L^k(x^*-x_0)| \]

因為\(\displaystyle \lim_{k \to \infty}L^k=0\)，故\(\displaystyle \lim_{k \to \infty}|x^*-x_{k+1}|=0\)(絕對值必定非負)，得\(x^*=x_{k+1}(k\to\infty)\)

3:迭代法的誤差式：

設某一次迭代後誤差為\(\epsilon\)，則可以推出：

\[|x_{k+1}-x_{k}| = |(x^*-x_k)-(x^*-x_{k+1})\geq|x^*-x_k|-|x^*-x_{k+1}|\geq|x^*-x_k|-L|x^*-x_k| \]

其中從\(|x^*-x_{k+1}| \leq L|x^*-x_k|\)則是因為對於一個可以收斂的函式而言，每一次迭代後，誤差總是減少的。

故可以輕鬆的計算出誤差估計式為：

\[|x^*-x_k|\leq\frac{1}{1-L}|x_{k+1}-x_k| \]

通過中值定理，又可以推出：

\[|x_{k+1}-x_k| = |\varphi'(\xi)(x_k-x_{k-1})| \]

因為\(|\varphi'(\xi)|\leq L\)，將上式遞推後放縮成第二種誤差估計式：

\[|x^*-x_k|\leq\frac{L^k}{1-L}|x_1-x_0| \]

其中L可以近似的看作是一個足夠接近方程解的函式導數的值。

3.3.3不動點迭代法的收斂速度

首先我們先直接給出一個結論：在不動點迭代過程中，第i步迭代和第i+1步迭代的表示式分別為\(e_i,e_{i+1}\)，且有\(\displaystyle \lim_{ i \to \infty}\frac{e_{i+1}}{e_i}=S<1\neq0\)，因此不動點迭代法是一個線性收斂的演算法。

如何能證明我們的演算法是一個線性收斂的函式呢？實際上非常簡單。證明過程如下：

\[若\varphi(x)連續可微,不動點的一個足夠近似的值為x_0，x_{i+1}和x_i是兩次迭代的估計值 \]

\[根據中值定理，\exists r\in (x_i,r),\frac{\varphi(x_i)-x_0}{x_i-x_0}=\varphi'(r) \]

\[\varphi(x_i) = x_{i+1},故x_{i+1}-r = \varphi'(r)(x_i-r) \]

\[e_{i+1}=\left|\varphi'(r)\right|e_i，即S=\left|\varphi'(r)\right|為常數 \]

我們可以在這裡舉幾個例子，例如我們希望求函式\(f(x) = x^3+x-1=0\)的根，倘若我們取\(\varphi(x) = 1-x^3\)，根據我們誤差分析的式子，包括我們收斂分析的公式來說，後一步的誤差將會是\(e_{i+1}=-3r^2e_i\)，我們預估的根大概是0.6823左右，可以發現我們的收斂速度的絕對值已經大於1，意味著每一步的誤差都在擴大，因此函式不會收斂。如果我們取\(\varphi(x) = \sqrt{1-x}\)，我們繼續計算出的收斂速度大約是0.7左右，很明顯誤差在以線性縮小。如果再次修改我們的遞推式\(\varphi(x) = \frac{1+3x^3}{1+3x^2}\)，我們的S幾乎是等於0的，因此也解釋了這種方法為何是最快收斂到我們的根。

//todo 圖片

再舉一個例子，函式\(f(x) = cosx-sinx=0\)的迭代求根過程，我先給定我們的近似值大約是0.7853982。這裡似乎很難構造我們的迭代遞推式，事實上只需要方程兩邊同時加x就可以湊出\(x=cosx-sinx+x=\varphi(x)\)。我們利用計算機編寫好程式並計算，最後得到了這樣一組迭代過程中的一些資料：

不難看出我們的\(\frac{e_i}{e_{i-1}}\)在一個恆定值0.414保持了相當的次數，這也是因為\(S=\left|\varphi'(r)\right|=|1-sinr-cosr|\approx0.414\)，印證了我們的推論是完全正確的，當我們不動點迭代在一個區間內，往往會以線性模式進行誤差的縮小。

迭代的終止

3.4牛頓法

牛頓在整個微積分、數值分析中都有著舉足輕重的地位，這裡闡述的牛頓法，全名牛頓-拉夫遜法，就是Taylor展開式的一部分。牛頓迭代法的核心思想就是：設法將一個非線性方程\(f(x)=0\)轉化為某種線性方程去求解。在數學意義上，我們知道泰勒公式可以將任意函式以簡單冪函式的形式表示出來，而在幾何意義上，我們是利用切線與X軸的交點去進行迭代，在根處，切線必過零點。若設\(f(x)=0\)的近似解為\(x_k\)，則方程可以用一階Taylor展開進行近似表示，牛頓迭代法的核心公式如下：

\[p_1(x) = f(x_k)+f^{'}(x_k)(x-x_k) \]

3.4.1普通牛頓法

從上述公式中我們知道，其中\(p_1(x)\)就是泰勒多項式的表達，將其看成是方程\(f(x)=0\)，通過迭代的思想，從而得到了一個線性方程：

\[x_{k+1} = x_k - \frac{f(x_k)}{f^{'}(x_k)} \]

這就是普通牛頓法的迭代公式。在幾何意義上，他就是一個切線與x軸的交點去逼近零點，如圖：

我們不斷的通過迭代這個過程，就能讓我們的值越來越精確。

3.4.2牛頓下山法

所有的迭代法都有一個無法逃過的缺點，如果選取的初始值離根太遠，則會導致迭代過程次數過多，並且有可能導致迭代發散，因為迭代法只具有區域性收斂性。

為了避免迭代失敗或時間過長，我們加上這樣一個條件用於保證數值穩定下降

\[|f(x_{k+1})|<|f(x_k)| \]

將這個條件和牛頓法結合再一起，即再保證數值穩定下降的過程中，利用牛頓法加快迭代，這就是牛頓下山法。具體的操作如下：

將牛頓法的結果\(\bar{x}_{k+1}\)與前一步的迭代值\(x_k\)進行加權平均作為新的迭代值\(x_{k+1}\)，

則有：

\[x_{k+1} = \lambda\bar{x}_{k+1} + (1-\lambda)x_k \]

或

\[x_{k+1} = x_k - \lambda\frac{f(x_k)}{f^{'}(x_k)} \]

其中\(\lambda(0\leq\lambda<1)\)稱為下山因子，它的值是一個逐步試探的過程，可以從1開始取值，一旦滿足\(|f(x_{k+1})|<|f(x_k)|\)則稱為下山成功，否則需要另選初始值\(x_0\)進行試算。

3.4.3簡單牛頓法

牛頓法可以說是一個非常有效的計算方法，準確度和迭代次數上都要比普通的迭代法要好得多，但是牛頓法最大的問題是我們需要求方程的導數，對於某些極其複雜的函式而言，導數是無法通過人工的方式計算，假如我們使用微積分的方式去求解導數，這對整個程式的效能會有比較大的影響，因此我們可以利用一個常數值\(\lambda\)來代替方程中的導數項。此時迭代公式為：

\[x_{k+1} =x_k - \frac{f(x_k)}{\lambda} \]

不過對於常數\(\lambda\)的取值是有限制的，因為我們需要保證迭代函式的收斂性，如果函式不收斂於\(x^*\)，那麼一切都沒有意義。於是有：

\[\varphi(x) = x - \frac{f(x)}{\lambda}\Longrightarrow\varphi^{'}(x) = 1-\frac{f^{'}(x)}{\lambda} \]

牛頓迭代法的收斂性遵循前文中普通迭代法的收斂性，於是可以得到：

\[|\varphi^{'}(x)| = |1-\frac{f^{'}(x)}{\lambda}|\Longrightarrow0<\frac{f^{'}(x)}{c}<2 \]

這樣我們就可以很輕鬆的確定下c的值了。

簡單牛頓法的幾何意義就簡單許多了，和我們之前討論的普通迭代法一致，只不過普通迭代法是將函式值和\(y=x\)進行處理變換，而簡單牛頓法則是和\(y= \lambda(x-x_k)+f(x_k)\)進行變換，本質是一致的。

///這裡只寫普通牛頓法，另外的函式由讀者自己思考
// 其中f1x為導數
public static double Newtown(Func<double, double> fx, Func<double, double> f1x, double x)
{
    var temp = f1x.Invoke(x);
    if (temp == 0)
    {
        throw new ArgumentException();
    }
    x = x - fx.Invoke(x) / temp;
    if (Math.Abs(fx.Invoke(x)) <= epsilon)
    {
        return x;
    }
    return Newtown(fx, f1x, x);
}

3.4.4 改進牛頓法——弦截法

弦截法是牛頓法的一種改進，保留了牛頓法中收斂速度快的優點，還克服了再牛頓法中需要計算函式導數值\(f^{'}\)的缺點。

弦截法中，我們用差商

\[\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} \]

去代替牛頓法中的導數值，於是可以得到以下離散化的迭代

\[x_{k+1} = x_k - \frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1}) \]

這種方法叫做雙點弦截法，如圖所示：

從圖中可以知道，弦截法一直利用兩點之間的連線作為迭代的內容，那麼，他的合理性在哪裡呢？

整個迭代法都離不開中值定理，這裡也是這樣，事實上，這個差商之所以可以對導數值進行替代，是因為中值定理中說過，連續函式中兩函式上的點的連線的斜率必為兩點之間某一點的導數值，並且迭代過程中，這兩點的中值會越來越逼近函式零點，事實上這已經說明了弦截法是牛頓法的改進方法了。

不過，如果將函式中\(x_{k-1}\)用一個定點\(x_0\)代替，這種方法叫做單點弦截法，幾何意義如圖所示：

//這裡就對雙點弦截法進行程式設計
public static double StringCut(Func<double, double> func, double x1, double x2)
{
    var temp = x1 - (func.Invoke(x1) / (func.Invoke(x1) - func.Invoke(x2))) * (x1 - x2);
    x2 = x1;
    x1 = temp;
    if (Math.Abs(func.Invoke(x1)) <= epsilon)
    {
        return x1;
    }
    return StringCut(func, x1, x2);
}

3.4.5 牛頓法的收斂性

牛頓法的收斂性有可能是一階收斂也有可能是二階收斂，具體要看我們的迭代函式。我們首先講講牛頓法二次收斂的情形。

假定我們牛頓法的迭代函式為\(g(x) = x-\frac{f(x)}{f'(x)}\)，那麼我們迭代函式的導數為\(g'(x) = \frac{f(x)f''(x)}{f'(x)^2}\)，設根為r，那麼易知\(g'(r)=0\)，迭代函式是收斂的。

我們再將\(f(r)\)進行泰勒展開，取\(x_i\)是一個足夠接近根的一次迭代結果，並且代入\(f(r)=0\)則有

\[f(r) = f(x_i)+(r-x_i)f'(x_i)+\frac{(r-x_i)^2}{2}f''(\xi),\xi \in (x_i,r) \]

\[f(r) = 0 =>x_i-\frac{f(x_i)}{f'(x_i)} -r = \frac{(r-x_i)^2}{2}\frac{f''(\xi)}{f'(x_i)} \]

若\(f'(x_i)\neq0\)對此函式迭代式在進行一次迭代可以得到以下式子

\[x_{i+1}-r=e_i^2\frac{f''(\xi)}{2f'(x_i)} \]

\[e_{i+1} = e_i^2\left|\frac{f''(\xi)}{2f'(x_i)}\right| \]

這樣我們就得到了牛頓法的二階收斂的定義，對於二階收斂，我們要求是需要在該處的一階導數值不為0。

不過我們的牛頓法並不是總為二次收斂的，在面臨多重根的時候，牛頓法往往會倒退到線性收斂的速度。，例如函式\(f(x)=x^m\)利用牛頓法求根的時候，牛頓公式如下：

\[x_{i+1} = x_i-\frac{x_i^m}{mx_i^{m-1}}=\frac{m-1}{m}x_i \]

唯一的根是0，但是0是一個m-1重根，因此得到

\[e_{i+1} = Se_i,S=\frac{m-1}{m} \]

我們這裡可以給出一個定義，若m+1階連續可微函式有一個m重根，那麼利用牛頓法會線性區域性收斂，誤差式如上。

收斂迭代的加速與精度

利用迭代法進行收斂的計算時，往往面臨一個很麻煩的問題，就是當執行足夠多次迭代的時候，我們的精度往往會不降反升，這是因為一種類似搖擺的問題出現，例如我規定此刻你站在座標1上，我希望你走到座標0上，你每次可以向左或者向右走三步，那麼第一次你向左走三步之後，得到的迭代結果時-2，反而離精確值更遠了。這是一個令人糟心的問題。這裡我們對這種現象進行一次解釋。

收斂的精度

y還是x？

事實上我們的誤差分為了前向誤差和後向誤差，什麼是前向誤差呢？實際上就是我們要求解的根與實際的差值，假設\(f(r)=0\)這個方程求根，前向誤差定義為\(\left|r-x_i\right|\)，而後向誤差指代的就是我們在函式值上的接近程度，也就是y軸的接近程度，我們定義為\(\left|f(x_i)|\)。

對於計算機而言，儲存的精度是有極限的，假定我們現在有一臺古老的計算機，他的儲存精度最多就到小數點後4位，那麼假設我們的迭代函式是\(f(x)=x^2\)，我們知道在精度範圍內，0.001是他能取到最小具有有效數字的精度解，如果我們現在通過迭代取到的數字是0.001的話，我們會發現一個驚人的問題，再進行一次迭代之後，得到\(f(0.001)=0.000001\)，由於我們的計算機無法儲存足夠的精度，也就是說這個值會被計算機看成是0，但我們知道這並不是我們的最優解。

這其實就是我們後向誤差足夠小，但是前向誤差並不夠小導致的。我們在實際的程式碼編寫的時候，一定要注意的一個地方就是我們的誤差在什麼時候應當利用前向誤差進行終止，什麼時候又應該利用後向誤差進行求解。

敏感性

讀者可以嘗試將迭代初始值設定成我們方程實際的根進行迭代，按著我們正常人的理解，他應該立即停止迭代直接返回我們的初始值即可。但是演算法往往是不盡人意的，在絕大多數工具或語言的內建求根函式中，即使我們傳入的初始值是根，往往得到的結果反而是一個不精確的解。例如著名的威爾金森多項式：

\[W(x) = (x-1)(x-2)...(x-20) \]

如果我們利用這個已經因式分解好了的式子進行求根運算，這似乎沒有必要，我們可以得到一個非常工整的方程根，他就是1-20之間所有的自然陣列成，一旦我們將其展開後投入運算，由於存在許多大數字，很多大數字會因為儲存、計算等過程中產生損失和誤差，這些誤差實際上是很小的數字，但是卻對我們的計算結果產生了很大的誤差，最終我們難以得到一個精確解。這就是對於求根過程中對數字擾動的敏感性問題。

我們進一步去探索一下究竟是什麼導致了誤差放大，並且能粗略的計算出這些誤差究竟有多少，我們假設是尋找一個函式\(f(x)=0\)的根r，同時對我們的函式做出一個很小的擾動\(\epsilon g(x)\)，其中我們的係數\(\epsilon\)很小，用\(\Delta r\)表示根的變化。那麼會有以下式子。

\[f(r+\Delta r)+\epsilon g(r+\Delta r) = 0 \]

將上式統統泰勒展開

\[f(r)+(\Delta r)f'(r)+\epsilon g(r)+\epsilon(\Delta r)g'(r)+O((\Delta r)^2) = 0,其中\displaystyle O((\Delta r)^2) \to 0,\Delta r \to 0 \]

忽略後面的高階無窮小，同時r又是根，於是可以得到一個關於r的變化值的近似

\[\Delta r (f'(r)+\epsilon g'(r)) \approx -f(r)-\epsilon g(r) =-\epsilon g(r) \]

又因為\(\epsilon\)也是一個很小的值（需要遠小於\(f'(r)\)），我們直接認為他就是0，得到了\(\Delta r\)的最終估計:

\[\Delta r \approx -\epsilon\frac{g(r)}{f'(r)} \]

這就是我們根的敏感性公式，利用這個公式，我們可以得出，如果出現了某些數字的擾動，根會如何變化。

這裡舉一個黃皮書上的例子，如果函式\(P(x) = (x-1)(x-2)...(x-6)-10^{-6}x^7\)，求出它的最大根。我們肉眼可見最大根應該是在6附近，如果沒有最後那一項，我們的最大根就是6，那麼問題就在於如果我去除最後一項，對根會有多大的影響呢？

那麼我們直接設近似函式為\(f(x)=(x-1)(x-2)...(x-6)\)，令我們的\(\epsilon g(x)=-10^{-6}x^7\)。通過我們的根敏感性公式計算一番之後，得到了

\[\Delta r \approx -\epsilon\frac{g(r)}{f'(r)} = 2332.8\times 10^{-6} \]

因此我們估計的最大根應當是6.0023328左右，讀者可以自行使用迭代法去求解最大根，最後得出的結果與這個估計值是極度接近的。這裡我們定義一個新東西誤差放大因子，也就是我們的\(\left|\frac{g(r)}{rf'(r)}\right|\)。

普通迭代加速

對於迭代過程，利用中值定理有

\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]

當\(\Delta x \rightarrow0\)，我們將\(\varphi^{'}(\xi)\)看成定值\(\lambda(\lambda<1)\)

於是有

\[\lambda(x^*-x_k) = x^*-x_{k+1}\longrightarrow x^* = \frac{1}{1-\lambda}\bar{x}_{k+1}-\frac{\lambda}{1-\lambda}x_k \]

故可以推出最終的迭代公式為

\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ x_{k+1} = \bar{x}_{k+1} + \frac{\lambda}{1-\lambda}(\bar{x}_{k+1}-x_k) \end{cases} \]

利用這種方式的好處就是再求得\(x_k\)時，利用加權的方式，使得迭代法變得有點類似像牛頓迭代法一樣變成切線性質的線而不是x=c或y=c的線，你經過畫圖和簡單的代數運算後，你會發現\(\bar{x}_{k+1}<x_{k+1}\)，也就是說達到了我們的加速目的。

Atiken加速法

Atiken加速法其實就是再普通加速迭代公式上在進行一次迭代，這裡直接寫出公式：

\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ \bar{\bar{x}}_{k+1} = \varphi(\bar{x}_{k+1})\\ x_{k+1} = \bar{\bar{x}}_{k+1} - \frac{(\bar{\bar{x}}_{k+1}-\bar{x}_{k+1})}{\bar{\bar{x}}_{k+1} - 2\bar{x}_{k+1}+x_k} \end{cases} \]

練習題

1. 試著不借助計算器計算\(\sqrt{2}\)的值
2. 試著證明\(f(x)=ax+b\)在牛頓法中會一步收斂
3. 程式設計題：由一個高度為10m的圓柱構成的發射井頂部是一個半球，體積\(400m^3\)，確定發射井底部班級
4. 設函式為\(f(x)=x^n-ax^{n-1},g(x)=x^n\)，利用敏感性公式手動估計\(f_\epsilon(x) = x^n-ax^{n-1}+\epsilon x^n\)的非零根

Reference

《數值分析》（原書第二版） —— Timothy Sauer

《數值計算方法》（清華大學出版社） —— 呂同富等

About Me

作　　者：WarrenRyan
出　　處：https://www.cnblogs.com/WarrenRyan/
本文對應視訊：BiliBili(待重錄)
關於作者：熱愛數學、熱愛機器學習，喜歡彈鋼琴的不知名小菜雞。
版權宣告：本文版權歸作者所有，歡迎轉載，但未經作者同意必須保留此段宣告，且在文章頁面明顯位置給出原文連結。若需商用，則必須聯絡作者獲得授權。
特此宣告：所有評論和私信都會在第一時間回覆。也歡迎園子的大大們指正錯誤，共同進步。或者直接私信我
聲援博主：如果您覺得文章對您有幫助，可以點選文章右下角【推薦】一下。您的鼓勵是作者堅持原創和持續寫作的最大動力！


博主一些其他平臺：
微信公眾號：寤言不寐
BiBili——小陳的學習記錄
Github——StevenEco
BiBili——記錄學習的小陳（計算機考研紀實）
掘金——小陳的學習記錄
知乎——小陳的學習記錄

聯絡方式

電子郵件：cxtionch@live.com

社交媒體聯絡二維碼：