十行代码实现牛顿方法

465 查看

从今天起,我将推出系列博文:十行代码解决数学问题。这第一篇我将介绍牛顿迭代法,这种古老的数值逼近法能够用来求解复杂多项式和任意可导函数的根。

首先假设我们有如下复杂多项式:

latex

现在我们想求得此多项式的根。然而根据伽罗瓦理论:对于五次及以上的多项式,没有统一的求根公式;因此我们不得不求助于数值方法。什么意思呢?在中学的时候我们就学过怎么解二次方程,同样地,对于三次或四次方程,我们也能得到其公式解。但是,伽罗瓦已经证明当多项式的次数在五次及以上的时候,不存在仅仅利用常规代数运算(包括加、减、乘、除)和根运算(包括开平方根,立方根等)就能得到多项式解的公式。

在 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 的计算公式:

latex (1)polynomial

综上所述,我们很容易写出 f(x) = 0 的近似求根算法,使其满足任意误差 e。很明显,在上例中,选取不同的起始估值点,可能找到不同的根。

现在你已经了解(牛顿迭代法)了,实现该算法的代码不超过10行(包括一空行代码),要是去掉打印结果的语句的话就更少了。为了使代码能够运行,我们还需要两个函数 f(x), f'(x),对于这个多项式的话,这两函数为:

现在我们就可以很容易的找到多项式的三个根了:

以上所有的代码,以及另一部分与 scipy.optimize.newton 方法比较测试的代码都可以在这个 Gist 上找到。还有别忘了,要是很难得到曲线的切线函数的话,可以使用 Sympy,代码在这里

虽然牛顿迭代法很强大,但在加速收敛过程的同时,可能会有一些问题,如果初始估值点选的不好的话甚至会导致其不收敛,比如这样。不管怎样,我希望你们能认识到牛顿迭代相当有用…大家有什么想法的话就发表在下面评论区吧。