上次讲解了如何使用 Mathematica 的 FindRoot[] 函数,数值求解代数方程、(常、偏)微分方程。不管是初值还是边值甚至是刚性问题,都可以使用该方法。本博客的目的不仅仅是介绍 Mathematica 中的函数,而是一步步剖析其中的原理,理解其中的算法。算法不仅可以用 Mathematica 来实现,以后有机会也将更新使用 Python、C++ 来实现我所介绍的所有算法,希望是订阅数超过 10000 之后。


FindRoot[]

我们看一下 Mathematica 的 FindRoot[] 函数使用了哪些算法,Method -> 可以指定方法,请查看帮助文档 tutorial/UnconstrainedOptimizationMethodsForSolvingNonlinearEquations

Method 解释
“Newton” 使用精确雅可比矩阵或者是有限差分近似在一个局部线性模型基础上的求解步骤
“Secant” 不需要使用导数进行求解(类比牛顿法),通过利用已完成的 n 个步骤建立一个雅可比的割线近似;需要在每个维度上都有两个初始条件
“Brent” 在一个维度上保持对根的包围(类比二分法);需要两个初始条件来包围一个根
“AffineCovariantNewton” 优化牛顿法,使用最少数量的雅可比矩阵的计算

牛顿法

牛顿法(Newton’s method),也称作牛顿-拉弗森方法(Newton-Raphson method),以 Isaac Newton 和 Joseph Raphson 的名字命名。牛顿法是一种非常有效的在实数域和复数域上近似求解方程的方法,在工程、物理、计算机科学等众多领域有着广泛的应用。在实际应用中,牛顿法广泛用于求解非线性方程、优化问题和数值分析中的根查找等领域。


原理

牛顿法的基本思想是将函数 f(x) 在某个初始近似值 处展开成泰勒级数,并取其线性部分来近似表示原函数。这个线性化后的方程容易求解,其解可以作为新的近似值 。然后,用这个新的近似值再去展开原函数,并重复这个过程,直到达到所需的精度或者函数值的变化小于某个阈值。

具体来说,对于函数 f(x) ,在 处的泰勒展开式可以写作:

其中, 是 f(x) 的导数。由于我们只取泰勒级数的前两项,这被称为线性泰勒级数近似。使用这个线性近似,我们可以得到一个线性方程:

解这个线性方程可以得到新的近似值

这个新的近似值 会比之前的近似值更接近 f(x) 的零点。


Mathematica

  • 使用牛顿法数值解方程 的近似根,以 1 这个根为例。首先令 ,画出 f(x) 的图形:
In[]:=
1
2
3
f[x_] = x^2 - x;
Plot[f[x], {x, -0.5, 2}]
dfdx[x_] = D[f[x], x]

dfdx 为 f(x) 的导函数。下面以 x=0.6 作为初始值进行迭代:

In[]:=
1
2
3
4
5
6
7
8
9
10
11
x0 = 0.6; error = Abs[f[x0]];
n = 0;
(*n 为迭代次数,为防止遇到不收敛的情况,将最大迭代次数设为 1000。
计算精度要求 0.0001,可以自行设定为需求精度*)
While[error > 0.0001 && n < 1000,
x0 = x0 - f[x0]/dfdx[x0];
(*Print[x0];*)
error = Abs[f[x0]];
n++]
n
x0
Out[]:=

可以看出经过5步迭代就收敛达到精度。


几何意义

牛顿法求解零点的过程,在几何上可以理解为寻找函数图像与 x 轴交点的横坐标。选取一个初始似值 ,便可以在函数图像上找到这一点,然后作切线,这条切线与 x 轴的交点就是下一个近似值 。不断重复这个过程,直到找到足够接近的零点。因此牛顿法也叫切线法

xf(x)-1-0.500.511.520.511.52f(x)=x2¡xx0x1x2x3


注意事项

  • 牛顿法在开始迭代之前需要取一个好的初始近似值。上面的例子如果选择 0.5 作为初始值就会遇到 1/0 这种情形。
  • 如果函数在零点附近是单调的,那么牛顿法会很快收敛。
  • 如果函数在零点附近有多个极值点,选择不当的初始近似值可能会导致迭代过程发散或不收敛。
  • 在实际应用中,可能需要对函数进行适当的修改或使用 modified Newton’s method 来改善收敛性或避免除以零的问题。

上面的例子仅包含了一个方程,下次会介绍如何利用牛顿法对多个方程求解。

请大家多多关注、点赞、收藏、转发,让更多有需要的人可以看到,以后会分享更多实用的算法。

万分感谢!🙏