數值方法導論/方程求根

求根

目標

• 找到二次方程和三次方程的解
• 推匯出公式並遵循以下方法求解非線性方程的演算法：
  • 二分法
  • 牛頓-拉夫森法
  • 割線法
  • 試位法

資源

• 關於各種方法的書籍章節
• 求根演算法
• 二分法
• 牛頓法
• 割線法
• 試位法

函式 f(x) 的根（或零點） 是產生輸出為 0 的 x 值。根可以是實數或複數。找到 $f(x)-g(x)$ 的根與解方程 $f(x)=g(x)$ 是相同的。解方程 是找到滿足方程指定條件的值。

低次（二次、三次和四次）多項式有閉式解，但數值方法可能更容易使用。為了解二次方程，我們可以使用二次公式

$ax^{2}+bx+c=0$

$x={\frac {-b\pm {\sqrt {b^{2}-4ac\ }}}{2a}}$

有許多 求根演算法 用於數值求解方程。

二分法

二分法從兩個猜測開始，並使用二分查詢演算法來改進答案。如果函式在兩個初始猜測之間是連續的，則保證二分法會收斂。這是一張說明該想法的圖片

二分法的優點包括在連續函式上保證收斂，並且誤差是有界的。

二分法的缺點包括收斂速度相對較慢，並且在某些函式上不收斂。

牛頓-拉夫森方法

牛頓法（又稱牛頓-拉夫森方法）是一種求解非線性方程的開放方法。與區間收縮法（例如二分法）相反，牛頓法只需要一個初始猜測，但它不保證收斂。

牛頓法的基本思想如下

給定一個“x”的函式 f 和一個初始猜測 $x_{0}$ 用於該函式的根，一個更好的猜測 $x_{1}$ 是

$x_{1}=x_{0}-{\frac {f(x_{0})}{f'(x_{0})}}$

這就是牛頓-拉夫森公式。

以下動畫說明了該方法

讓我們使用牛頓法求解以下方程

$f(x)=x^{2}+2$

$f(x)=sin(x)=0$

$f(x)=a-{\frac {1}{x}}=0$

from sympy import Symbol, diff, sin
x = Symbol('x')
y = sin(x)
derivative = diff(y, x)

'''
Solve f with initial guess g using Newton's method.
Stop when the absolute relative approximate error is
smaller than the tolerance or the max # iterations
is reached.
'''
def newton(f, derivative, g, tolerance, max_iteration):
  x_previous = g
  for i in range(max_iteration):
    x_current = x_previous - \
      y.subs(x, x_previous).evalf()/derivative.subs(x, x_previous).evalf()

    error = abs((x_current-x_previous)/x_current)
    print "current x:", x_current, " error:", error

    x_previous = x_current

    if error < tolerance:
      break;

newton(y, derivative, 5, 0.005, 15)

$ python newton.py 
current x: 8.38051500624659  error: 0.403377955140806
current x: 10.1008867567293  error: 0.170318883075941
current x: 9.29864101772707  error: 0.0862755898924158
current x: 9.42545121429349  error: 0.0134540186653470
current x: 9.42477796066766  error: 7.14344283379346e-5

牛頓法的侷限性

除以零

由於牛頓公式中涉及除法，當分母為零時，該方法將無法繼續正確執行。以下程式演示了此問題。

from sympy import Symbol, diff
x = Symbol('x')
y = 4 - 1.0/x
derivative = diff(y, x)

'''
Solve f with initial guess g using Newton's method.
Stop when the absolute relative approximate error is
smaller than the tolerance or the max # iterations
is reached.
'''      
def newton(f, derivative, g, tolerance, max_iteration):
  x_previous = g
  for i in range(max_iteration):
    x_current = x_previous - \
      y.subs(x, x_previous)/derivative.subs(x, x_previous)

    error = abs((x_current-x_previous)/x_current) 
    print "current x:", x_current, " error:", error
    
    x_previous = x_current

    if error < tolerance:
      break;

newton(y, derivative, 0.5, 0.005, 15)

輸出

$ python newton.py 
current x: 0  error: oo
current x: nan  error: nan
current x: nan  error: nan
current x: nan  error: nan
current x: nan  error: nan
current x: nan  error: nan
...

拐點處的發散

當初始猜測接近拐點時，該方法可能會偏離所需的根。

以下程式碼使用牛頓法的相同實現（假設newton函式已匯入）演示了拐點（x=0）處的發散。

from sympy import Symbol, diff
x = Symbol('x')
y = (x-1)**3 + 0.5

derivative = diff(y, x)
newton(y, derivative, 5, 0.005, 15)

$ python newton.py 
current x: 3.65625000000000  error: 0.367521367521368
current x: 2.74721164936563  error: 0.330894909696631
current x: 2.11021215705302  error: 0.301865141940137
current x: 1.60492272680972  error: 0.314837232847788
current x: 0.947823175470885  error: 0.693272298403532
current x: -60.2548036313237  error: 1.01573025084058
current x: -39.8365801731819  error: 0.512549605648313
current x: -26.2244867245776  error: 0.519060433539279
current x: -17.1498826852607  error: 0.529135050417335
current x: -11.1004277326071  error: 0.544974941360464
current x: -7.06809009701998  error: 0.570498901434099
current x: -4.38128712811798  error: 0.613245124168828
current x: -2.59328016409484  error: 0.689476975445585
current x: -1.40842833626215  error: 0.841258157995587
current x: -0.634351912031382  error: 1.22026340513730

割線法

割線法類似於牛頓法，因為它是一種開方法，並使用交點來獲得根的改進估計。割線法透過使用割線的斜率來估計導數值，從而避免計算一階導數。

下圖說明了割線法，其中有兩個初始猜測和兩個連續的改進估計：

割線法的公式如下

$x_{i+1}=x_{i}-{\frac {f(x_{i})(x_{i}-x_{i-1})}{f(x_{i})-f(x_{i-1})}}$

試位法

試位法類似於二分法，因為它需要兩個初始猜測（區間法）。試位法不是使用中點作為改進的猜測，而是使用穿過兩個端點的割線的根。下圖說明了這個想法。

試位法的公式如下

$x_{r}={\frac {x_{u}f(x_{l})-x_{l}f(x_{u})}{f(x_{l})-f(x_{u})}}$