十行程式碼實現牛頓方法

從今天起，我將推出系列博文：十行程式碼解決數學問題。這第一篇我將介紹牛頓迭代法，這種古老的數值逼近法能夠用來求解複雜多項式和任意可導函式的根。

首先假設我們有如下複雜多項式：

現在我們想求得此多項式的根。然而根據伽羅瓦理論：對於五次及以上的多項式，沒有統一的求根公式；因此我們不得不求助於數值方法。什麼意思呢？在中學的時候我們就學過怎麼解二次方程，同樣地，對於三次或四次方程，我們也能得到其公式解。但是，伽羅瓦已經證明當多項式的次數在五次及以上的時候，不存在僅僅利用常規代數運算（包括加、減、乘、除）和根運算（包括開平方根，立方根等）就能得到多項式解的公式。

在 reddit 上有好心的朋友說我沒有真正理解伽羅瓦理論。很抱歉，畢竟我不是一位數學家。實際上，我們可以將上面的多項式因式分解成 x^2(x-1)(6*x^2 + x – 3)，然後再用傳統的三次求解法解出。此外，我所引用的是阿貝爾-魯菲尼定理(Abel-Ruffini Theorem)，該定理僅用來求解五次及以上的一般多項式。儘管如此，這個例子還是有用的，我們能夠用它來說明如何應用牛頓迭代法，而對於其他的多項式，則依此類推。

在這種情況下，我們必須用線性逼近的方法來解決。最簡單的就是介值定理：假設函式 f(x) 在閉區間 [a, b] 上連續，且有 f(a) < y < f(b)，那麼在 [a, b] 內至少存在一點 x，滿足 f(x) = y。

因此根據介值定理，我們可以找到兩點 x1、x2 ，滿足 f(x1) > 0 且 f(x2) < 0，那麼使得 f(x) = 0 的點必定在 x1 和 x2 之間。接下來，我們再觀察 x3 = (x2 – x1)/2，如果 f(x3) > 0，那麼令 x4 = (x2 – x3)/2，再判斷 f(x4) 是否小於0；依次類推進行下去，可以得到 xn，使得 f(xn) 越來越趨近於 0。然而這種方法收斂很慢，故牛頓設計了一種更好的方法來加速收斂過程（假設能找到根的話）。

該多項式的曲線如下圖所示，可以看出，該多項式有三個根，分別在 x 等於 0、1，以及 0 到 1 之間的位置。那麼我們怎樣才能得到這些根呢？

在牛頓迭代法中，先任選一點 x0 作為根的近似估值，過點 (x0, f(x0)) 做曲線的切線，那麼切線的斜率即為 f'(x0)。然後取切線與 x 軸的交點 x1 作為新的估值點，重複上一步驟可以計算得到 f'(x1)、x2。以此類推，我們可以找到一點 xn，使得 |f(xn)| 趨於0，也就是說如果我們將上述操作進行下去的話，那麼我們離多項式的根就越來越近了。因此，我們只要得到了 xn 點處曲線的切線斜率 f'(xn)，再由此切線與 x 軸的交點便可以匯出 x n+1 的計算公式：

綜上所述，我們很容易寫出 f(x) = 0 的近似求根演算法，使其滿足任意誤差 e。很明顯，在上例中，選取不同的起始估值點，可能找到不同的根。

現在你已經瞭解(牛頓迭代法)了，實現該演算法的程式碼不超過10行(包括一空行程式碼)，要是去掉列印結果的語句的話就更少了。為了使程式碼能夠執行，我們還需要兩個函式 f(x), f'(x)，對於這個多項式的話，這兩函式為：

現在我們就可以很容易的找到多項式的三個根了：

以上所有的程式碼，以及另一部分與 scipy.optimize.newton 方法比較測試的程式碼都可以在這個 Gist 上找到。還有別忘了，要是很難得到曲線的切線函式的話，可以使用 Sympy，程式碼在這裡。

雖然牛頓迭代法很強大，但在加速收斂過程的同時，可能會有一些問題，如果初始估值點選的不好的話甚至會導致其不收斂，比如這樣。不管怎樣，我希望你們能認識到牛頓迭代相當有用…大家有什麼想法的話就發表在下面評論區吧。