数值方法入门

• 下载源代码 - 90.3 KB

引言

我最近将《数学函数教程：第2部分》中的递归下降解析器添加到了我的《多项式项目》中。在添加的同时，我还将其制作成了一个库，将Main()中的测试替换为名为UnitTestProjectMath的单元测试项目，并将整个项目更新为.NET 4.5.2。

不久之后，我在当地图书馆的书展上偶然发现了一本有趣的《数值数学与计算》（Ward Cheney和David Kincaid合著，Cengage Learning出版社，2012年，第7版）。最先吸引我眼球的是计算多项式的霍纳方法（第8页）。然而，书中只提供了伪代码，所以我认为我会写一些C#代码来实现该算法，而一个合适的地方是《多项式项目》，因为该项目已经处理了多项式。虽然Cheney和Kincaid的书籍涵盖了许多主题，但我将书中提供的伪代码转换成了C#，涵盖了以下六个主题，并通过以下链接进行了详细介绍：

1. 浮点数表示
2. 霍纳方法用于计算多项式
3. 二分法求根
4. 高斯消元法解方程组
5. 多项式插值
6. 近似计算圆周率

Using the Code

下载并解压项目，然后使用Visual Studio打开Polynomial.sln。它包含以下项目：

1. ComputePI - 一个Windows Forms项目，用于演示如何通过嵌入三角形来计算圆周率，如下图所示。
2. NumericalMath - 用于实现书中伪代码的C#代码
3. PolynomialLib - 最初的《多项式项目》已转换为库
4. RDPLib - 递归下降解析器 (RDP)。这与《数学函数教程：第2部分》中的RDP相同。
5. UnitTestProjectMath - 一个用于PolynomialLib和NumericalMath的单元测试项目

按F6键构建解决方案。

浮点数表示

以下是第一章中的一段引言，总结了理解浮点运算局限性的重要性：

人类历史上从未有过如此快速地产生如此多错误答案的可能性。.

任何学习过入门计算机科学课程的人，当他们运行一段类似这样的代码时，都会熟悉计算机上浮点数学的局限性。

	double sum = 0.0;
	for (int i = 0; i < 10; i++)
		sum += 0.1;

在这个例子中，sum的值不会精确地等于数学上预期的1.0。这是因为0.1并没有被精确存储。使用Jon Skeet的ToExactString(0.1)方法，我们可以看到0.1实际上存储为：

0.1000000000000000055511151231257827021181583404541015625

在NumericalMath项目中的FloatingPoint.cs文件中，我展示了如何存储双精度浮点数。例如，十进制数-3500.0以64位存储：一个符号位、11位的指数和52位的尾数。

1 10000001010 1011010110000000000000000000000000000000000000000000

通过将10000001010₂转换为其十进制等价值1034₁₀，然后从中减去1023₁₀，计算出指数，结果为11₁₀。尾数中的第一位（1）表示2^-1，我们忽略第二位（零），第三位表示2^-3，第四位表示2^-4，依此类推，将所有分数相加结果为0.708984375₁₀。正如在NumericalMath FloatingPoint.cs的ToBinaryString()中所展示的，将其加上1.0得到1.708984375₁₀。使用指数11₁₀，我们然后将2¹¹（即2048₁₀）乘以1.708984375₁₀，得到3500₁₀。由于符号位是1，最终结果是-3500₁₀。

霍纳方法用于计算多项式

Cheney和Kincaid在第8页描述了一种称为霍纳方法（或霍纳法则）的技术，用于仅使用加法和乘法（不进行指数计算）来计算多项式。代码很简单；我将提到我在NumericalMath HornersMethod.cs中提供了两个版本：HornerLINQ()和Horner()。后者（非LINQ版本）速度稍快，正如“测试资源管理器”窗口中的结果所示。

请参阅UnitTest1.cs中的TestHornerLINQ()和TestHorner()。

二分法求根

Cheney和Kincaid在第114页讨论了一种在给定区间内寻找连续函数根的方法。通过数学函数教程：第2部分，我们可以从下面的图中看到，方程f(x) = x^3.0 - 3.0 * x + 1.0在区间[0, 1]内的根约为0.34。

在BisectionMethod.cs的Bisection()方法中实现了Cheney和Kincaid的技术，以下是找到该区间内上述方程根的结果（有关函数f(x)，请参阅UnitTest1.cs中的MyFunction()；有关输出此表中数字的代码，请参阅NumericalMath BisectionMethod.cs中的#region debug）：

迭代	根	f(根)	Error（错误）
1	0.50000	-3.75e-001	5.00e-001
2	0.25000	2.66e-001	2.50e-001
3	0.37500	-7.23e-002	1.25e-001
4	0.31250	9.30e-002	6.25e-002
5	0.34375	9.37e-003	3.13e-002
6	0.35938	-3.17e-002	1.56e-002
7	0.35156	-1.12e-002	7.81e-003
8	0.34766	-9.49e-004	3.91e-003
9	0.34570	4.21e-003	1.95e-003
10	0.34668	1.63e-003	9.77e-004

由于Epsilon（最大绝对误差值）为0.001，经过10次迭代后找到的f(x)的根为0.34668。

请参阅UnitTest1.cs中的TestBisectionMethod()。

高斯消元法解方程组

在第71页，Cheney和Kincaid讨论了使用一种称为高斯消元法的技术来求解线性方程组。他们在第85页介绍了一种改进方法，称为部分主元法。GaussianElimination.cs中的GaussSolve()方法展示了这种技术。这是一个两阶段的过程：

1. 使输入矩阵进入上三角形式
2. 然后进行回代

上三角形式为（请参阅NumericalMath GaussianElimination.cs中的#region debug，其中包含输出此内容的的代码）：

12.00x₁ + -8.00x₂ + 6.00x₃ + 10.00x₄ = 26.00
0.00x₁ + -11.00x₂ + 7.50x₃ + 0.50x₄ = -25.50
0.00x₁ + 0.00x₂ + 4.00x₃ + -13.00x₄ = -21.00
0.00x₁ + 0.00x₂ + 0.00x₃ + 0.27x₄ = 0.27

现在，我们从底部开始回代。因此：

0.27 * x₄ = 0.27，x₄ = 1。

然后第三个方程变为：

4.00x₃ + -13.00 = -21.00，或者
4.00x₃ = +8.00，x₃ = +2.00，以此类推，计算出x₂和x₁，分别为-2和1。

请参阅UnitTest1.cs中的TestGaussianElimination()。

多项式插值

Cheney和Kincaid在第153页描述了一种多项式插值方法，并在NumericalMath PolynomialInterpolation.cs中实现了它。例如，给定一个如下所示的水温和粘度对照表，8°C时的粘度是多少？

温度 °C (x)	粘度 mPa·s (y)
0.0	1.792
5.0	1.519
8.0	?
10.0	1.308
15.0	1.140

首先，我们通过调用CalculateCoefficients()，传入表中x和y值以及用于接收系数的数组，来计算三次多项式的系数。然后将传递给CalculateCoefficients()的double[] coefficients数组，与x值和8（需要y值的x值）一起传递给Evaluate()。返回的y值为1.38。使用简单的线性插值……

double viscosity = 1.519 - ((8.0 - 5.0) * (1.519 - 1.308) / (10.0 - 5.0));

……结果为1.39，这让我们对结果有信心。请参阅UnitTest1.cs中的TestPolynomialInterpolation()。

近似计算圆周率

作为最后一个例子，我想做一些更有趣的事情。Cheney和Kincaid提供了一种古代数学家用来近似计算圆周率的方法（第35页），即在单位圆内嵌入三角形并计算它们的面积。其中一位这样做的数学家是阿基米德。在ComputePI项目中，我展示了一个Windows Forms项目来做这件事。当您选择迭代次数时，圆周率的近似值、绝对误差和相对误差（Cheney和Kincaid，第6页）会按如上图所示进行计算。我相信阿基米德本人会对自己能用3GHz四核处理器和一些代码做到的事情印象深刻。

结论

现代数字计算机能够比人们（甚至几十年前）想象的可能的速度更快地解决数学问题。然而，由于浮点数表示的局限性，正如开头引言中所明确指出的那样，一个人必须始终对结果保持谨慎，因为由于舍入误差、截断误差或迭代方法未能收敛到解，答案可能从非常准确到极其不准确。

历史

• 版本 2.0.0.0