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

在从二分逼近领略计算科学的魅力一文中，我们介绍了单调函数的求根公式（有零点），如 $f(x)=2x^2+3.2x-1.8$。我们能否采用二分逼近的原理，求解一个数的平方根（$x^2=n$）呢，自然地，我们将 $x^2=n$ 转换为求解 $f(x)=x^2-n$ 的零点，也即根，在 $[0,n]$ 的区间上，$f(x)=x^2-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]$ 之间采用二分逼近的方法查找，而 $\sqrt 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)$ 零点的 $x_0$，计算相应的 $f(x_0)$和切线斜率$f'(x_0)$（这里 $f'$ 表示函数 $f$ 的导数）。然后我们计算穿过点 $(x_0, f(x_0))$ 并且斜率为 $f'(x_0)$ 的直线和 $x$ 轴的交点的 $x$ 坐标，也就是求如下方程的解：

$$f(x_0)=(x_0-x)f'(x_0)$$

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

$$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$$

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

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

$$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}=x_n-\frac{f(x_n)}{2x_n}$$

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)-x^3=0$$
我们令 $f(x)=\cos(x)-x^3$，两边求导得， $f'(x)=-\sin(x)-3x^2$，由于 $-1\leq\cos(x)\leq1$（对于所有 $x$），则 $-1\leq x^3\leq 1$，即 $-1\leq x\leq 1$，可知方程的根位于 $[0,1]$，我们从 $x_0$ 开始：

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] 牛顿法