C# 中的三种求根方法

• 下载源代码 - 9.35 KB

引言

应用程序经常需要求解非线性方程。该问题的标准形式是求解一个数 x，使得 f(x) = 0。本文介绍了三种求解此问题的算法，并阐述了它们的优缺点。这些算法是用 C# 4.0 实现的。

背景

最著名的两种求根算法是二分法和牛顿法。简而言之，前者速度慢但稳定，后者速度快但不稳定。布伦特法既稳定，速度通常也比二分法快得多。

二分法完全可靠。假设你知道 f(a) 是负的，而 f(b) 是正的。只要 f 是一个连续函数，那么在 a 和 b 之间必然存在一个值 x，使得 f(x) = 0。二分法保证能在指定的容差范围内找到 x。事实上，在二分法的 n 步之后，误差 |f(x)| 保证小于 |f(b) - f(a)| 2^-n。粗略地说，该方法在每一步都能提供一个位的精度。

存在一类函数（例如，凸函数），对于这类函数，从任何初始猜测值开始，牛顿法都保证收敛。然而，在实践中，当没有理论保证其能正常工作时，该方法也经常被使用。这种方法通常会成功。但当牛顿法失败时，它可能以非常糟糕的方式失败。每次迭代不是接近期望的解，而是离解越来越远。如果被优化的函数近似于平坦，该方法可能一步就跳到离解非常远的地方。

当牛顿法起作用时，它效果非常好。一旦它开始收敛，结果中的正确位数会在每一步翻倍。而二分法在每一步都将误差减半，牛顿法则是将误差平方。

我曾见过一本教科书，它通过引用一首描绘牛顿法行为的童谣来介绍牛顿法。

从前有个女孩
她额头正中间
有个大卷发。
当她乖的时候
她非常非常乖，
但当她坏的时候，她就非常糟糕。.

布伦特法和二分法一样，保证收敛。在最坏的情况下，它可能比二分法稍慢。然而，它的收敛速度通常比二分法更接近牛顿法。如果函数是“良好”的（技术上讲，其导数是 Lipshitz 连续的），那么每一步的误差会在下一步中被提高到 1.6 次方。（确切的幂是黄金分割比例：(1 + v5)/2 = 1.1618...）

为了比较这三种方法的收敛速度，您可以像下面这样运行这三种方法。假设您为每种方法提供一个初始猜测值，您的初始猜测值精度在 1/2 以内，也就是说，一个位。

迭代一次后，二分法的答案将精确到两位。再迭代一次后，答案将精确到三位。位精度序列可能是 1, 2, 3, 4, ...

如果您的初始猜测值足够接近解，牛顿法的精度将是 1, 2, 4, 8, 16, ... 位。

布伦特法的精度可能是 1, 2, 3, 5, 8, 13, 22, ... 位。（您可能会注意到这些是斐波那契数。这并非巧合。连续斐波那契数的比率约等于黄金分割比例，并且随着数字的增大，这个近似值会提高。）

对于许多应用来说，布伦特法是理想的选择。它与二分法一样稳定，但速度却能与牛顿法媲美。布伦特法相对于牛顿法的另一个优点是，前者不需要导数函数作为参数，而后者需要。（布伦特法为了保证快速收敛，要求导数存在，但算法并不要求您编写一个函数来计算导数。）

Using the Code

RootFinding 类具有所有 static 方法；无需实例化对象即可使用其方法。该类有六种方法，每种方法（Bisect、Brent 和 Newton）都有两个重载版本。每个算法名称对应两个重载方法，一个简化版和一个详细版。详细版方法有两个额外的参数：iterationsUsed，即用于生成解的算法步数，以及 errorEstimate，即解的剩余误差的估计值。

求根的数学理论总是解决 f(x) = 0 形式的问题。如果要解一个像 x + sin(x) = 8 这样的方程，您首先需要定义一个新的函数 g(x) = x + sin(x) - 8，然后解 g(x) = 0。这种标准化在理论上是可行的，但在实践中有点麻烦，因为人们通常希望求解 f(x) = c 形式的方程，其中 c 是某个常数。为了方便起见，RootFinding 中的方法允许用户指定一个名为 target 的值，它代表这个常数 c，默认值为 0。

布伦特法的代码摘自该算法的创建者 Richard P. Brent 的著作《无导数最小化算法》(Algorithms for Minimization Without Derivatives)。这本书用 Algol 60 语言给出了该算法的源代码。这段代码被翻译成了 C#，并保留了原始的变量名。甚至还有一个 goto 语句。我选择保持与原始代码的一致性，而不是采用更现代的代码风格。

以下是使用布伦特算法简化接口的示例。

// Solve f(x) = 0 to default tolerance. 
// Look between -2 and 1.5. 

double root = RootFinding.Brent( new FunctionOfOneVariable(f), -2, 1.5);

以下是使用同一算法的更详细接口的示例。

// Solve f(x) = 0.2 to 10 decimal places.
// Look between x = 1 and x = 5.

double root = RootFinding.Brent(
    new FunctionOfOneVariable(f), // function to find root of, cast as delegate
    1.0,                          // left end of bracket
    5.0,                          // right end of bracket 
    1e-10,                        // tolerance 
    0.2,                          // target 
    out iterationsUsed,           // number of steps the algorithm used
    out errorEstimate             // estimate of the error in the result
);

关注点

布伦特法和二分法都需要起始值 a 和 b，使得函数 f(x) 在这两个点的值符号相反。如何找到这些起始点？您通常可以从问题的上下文中知道良好的起始点。如果您知道这些点一定存在，但不知道具体值，您可以搜索起始值。例如，您可以从 a = b 开始，然后逐渐减小 a 并增大 b，直到 f(a) 和 f(b) 符号相反。

历史

• 2010 年 5 月 6 日：初始发布