[數學]求函數近似零點的方法：二分法與牛頓法

問題

給定一個連續單變量函數\(f(x)\)，求這個函數的零點\(x_0\)。要求可控制誤差。

解決方案

二分法與牛頓法都是適合計算機的解決方案。不過，牛頓法遠快於二分法，寫起來也更簡單，但是更難理解。

二分法

算法是這樣的：

1. 找出（不管用什么方法，甚至看圖像也行）兩個值：\(l\)(low)與\(h\)(high)。使得\(f(l)<0\)且\(f(h)>0\)。由於連續性，\(x_0\)在\(l\)與\(h\)之間。
2. 設\(m=\frac{l+h}{2}\)，判斷\(f(m)\)的正負性：
   \(f(m)>0 \Leftrightarrow x_0在l與m之間\)
   \(f(m)=0 \Leftrightarrow x_0=m\)
   \(f(m)<0 \Leftrightarrow x_0在m與h之間\)

以上便是一次迭代。每經過一次迭代，\(x_0\)的范圍就會縮小一半。於是，經過多次迭代，\(l \mapsto h\)，我們就能將\(l\)或\(h\)作為\(x_0\)的近似值。

下面用二分法求函數\(f(x)=3^{x+1}+x^3\)的零點：

1. 輸入函數；畫圖
   
   由函數圖像可以猜測\(l=-10\)與\(h=0\)
2. 寫個程序進行迭代（誤差設定為\(10^{-5}\))
   
   花了二十次。

牛頓法

牛頓法的算法是這樣的：

假設\(x\)是猜測的\(x_0\)值，且\(x_1=x-\frac{f(x)}{f \prime (x)}\)，那么\(x_1\)比\(x\)更接近\(x_0\)。

以我現有的知識，無法理解牛頓法。但寫個程序還是可以的：

1. 我們已經畫好了圖，由圖像可以看到\(0\)比\(-10\)更接近\(x_0\)，所以我們猜測\(0 \mapsto x\)
2. 進行迭代（誤差一樣，設定為\(10^{-5}\))
   
   可以看到，牛頓法的代碼更精簡，而且更快：只需要5次。

結語

牛頓法效率遠高於二分法。但是牛頓法公式是怎么求出來的呢？希望在以后的學習中能夠了解。

最后，附上兩段Eigenmath代碼，是經過排版與添加注釋的二分法與牛頓法演示代碼：
已經在Eigenmath 137上測試通過。

二分法

f(x) = 3^(x+1)+x^3                       
draw(f)
l=-10                               # low
h=0                                 # high
cnt=0                               # count, 記錄迭代次數
error=10^(-5)                       # 誤差。若abs(high - low)<error則輸出結果
for(k,1,10^9,                       # 無限循環，直到小於誤差
    do(
        cnt=cnt+1,                  # 記錄迭代次數
        m=(l+h)/2,                  # 取二分點m
        test(
            f(m)>0,                 # 若f(m)>0
            h=m,                    #   h=m
            l=m                     # 否則 l=m
        ),
        test(
            abs(h-l) > error,       # 若結果大於誤差
            1,                      #   什么都不做，繼續迭代
            do(                     # 否則
                print(
                    float(l),       #   輸出low（小數形式）
                    float(h),       #   輸出high（小數形式）
                    cnt             #   輸出迭代次數
                ),
                stop                #   退出迭代
            )
        )
    )
)

牛頓法

f(x) = 3^(x+1)+x^3
draw(f)
guess=0                                                         # 猜測的x0的值
cnt=0                                                           # 記錄迭代次數（同上）
error=10^(-5)                                                   # 誤差（同上）
for(k,1,10^9,                                                   # 無限循環（同上）
    do(
        cnt=cnt+1,                                              # 記錄迭代次數（同上）
        x1=float(guess- f(guess)/eval(d(f,x),x,guess)),         # 使用牛頓法公式得出一個比guess更接近x0的值
        test( abs(x1-guess) > error,                            # 若猜測值增量大於誤差（表示還沒猜到）
            guess=x1,                                           #   更新猜測值guess，繼續迭代
            do(                                                 # 否則
                print(x1, cnt),                                 #   輸出接近x0的猜測值與迭代次數
                stop                                            #   退出迭代
            )
        )
    )
)