数值方法导论/方程求根

求根

目标

• 找到二次方程和三次方程的解
• 推导出公式并遵循以下方法求解非线性方程的算法：
  • 二分法
  • 牛顿-拉夫森法
  • 割线法
  • 试位法

资源

• 关于各种方法的书籍章节
• 求根算法
• 二分法
• 牛顿法
• 割线法
• 试位法

函数 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})}}$