平方根的计算（二分逼近、牛顿拉普生法）

在从二分逼近领略计算科学的魅力一文中，我们介绍了单调函数的求根公式（有零点），如 f(x)=2x2+3.2x−1.8。我们能否采用二分逼近的原理，求解一个数的平方根（x2=n）呢，自然地，我们将 x2=n 转换为求解 f(x)=x2−n 的零点，也即根，在 [0,n] 的区间上，f(x)=x2−n 自然单调递增，在此观点上我们几乎不需改变原始的求根函数，只需重新定义 f(x) 的形式即可。

typedef double (*FuncPtr)(double);double root(double x1, double x2, FuncPtr f){    double x = (x1+x2)/2;    while (abs(f(x)) > 1e-9)    {            // 1e-9：表示epsilon，一个小量        f(x) > 0 ? x2 = x : x1 = x;        x = (x1+x2)/2;    }    return x;}template<int N>double f1(double x){                    // 非类型模板参数仅限于常整数（包括 enum 枚举类型）                    //或者指向外部链接对象的指针    return x*x - N;}double f2(double x){    return 2*x*x+3.2*x-1.8}int main(int, char**){    std::cout << root(0, 5, f1<5>) << std::endl;                    // 5 的平方根    std::cout << root(-.8, .8, f2) << std::endl;                    // f(x)=2*x^2+3.2x-1.8 在区间 [-.8,.8] 内的零点    return 0;}

关于 C++ 模板编程的非类型模板参数，请参见 C++基础——非类型模板参数 。

是不是很完美，也不是，二分逼近求解实数的平方根时存在一个隐蔽的 bug，也即是要开平方的实数 x 小于1时，程序将永远无法退出，陷入死循环，因为我们是在 [0,x] 之间采用二分逼近的方法查找，而 x√>x，当 x<1 时，也即不在 [0,x] 的区间之内。一种精巧的改进即是，在调用端增加一个判断逻辑：

double n = .25;double f3(double x){    return x**x-n;}int main(int, char**){       std::cout << root(0, n>1 ? n : 1, f3) << std::endl;    return 0;}

牛顿拉普生方法

牛顿法（Newton’s method）又称为牛顿-拉弗森方法（Newton-Raphson method），它是一种在实数域和复数域上近似求解方程的方法。方法使用函数 f(x) 的泰勒级数的前面几项来寻找方程 f(x)=0 的根。该法效率极高，应用极广，并不限于求解实数的平方根，相反求解实数的平方根只是其一个具体的应用而已。

接下来我们介绍一个更具效率的平方根的计算方法，即为牛顿拉普生方法（Newton-Raphson method）。

方法说明：
首先，选择一个接近函数 f(x) 零点的 x0，计算相应的 f(x0)和切线斜率f′(x0)（这里 f′ 表示函数 f 的导数）。然后我们计算穿过点 (x0,f(x0)) 并且斜率为 f′(x0) 的直线和 x 轴的交点的 x 坐标，也就是求如下方程的解：

f(x0)=(x0−x)f′(x0)

我们将新求得的点的x坐标命名为 x1，通常 x1 会比 x0 更接近方程 f(x)=0 的解。因此我们现在可以利用 x1 开始下一轮迭代。迭代公式可化简为如下所示：

xn+1=xn−f(xn)f′(xn)

已经证明，如果 f′ 是连续的，并且待求的零点 x 是孤立的，那么在零点 x 周围存在一个区域，只要初始值 x0 位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果 f′(x) 不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

对于求解实数平方根的函数 f(x)=x2−n， 自然其根的迭代公式为：

xn+1=xn−f(xn)f′(xn)=xn−f(xn)2xn

double newton(double x, FuncPtr f){    while (abs(f(x))>1e-9)        x -= f(x)/(2*x);    return x;}double f(doubel x){    return x*x-.25;}int main(int, char**){    std::cout << newton(1., f) << std::endl;    return 0;}

牛顿拉普生方法的一个应用

求解如下方程的根：

cos(x)−x3=0

我们令 f(x)=cos(x)−x3，两边求导得，f′(x)=−sin(x)−3x2，由于 −1≤cos(x)≤1（对于所有 x），则 −1≤x3≤1，即 −1≤x≤1，可知方程的根位于 [0,1]，我们从 x0 开始：

double f(double x){    return cos(x)-x*x*x;}double f_prime(double x){    return -sin(x)-3*x*x;}double newton(double x, FuncPtr f){    while (abs(f(x)) > 1e-9)    {        x -= f(x)/f_prime(x);        std::cout << x << std::endl;    }    return x;}int main(int, char**){    std::cout << newton2(.5, f) << std::endl;    return 0;}

References

[1] 牛顿法