Python代码中的数学之美：用牛顿逼近法计算2的算术平方根

Posted 2021-05-08 天元浪子

编程的核心是算法，而算法的核心是数学——单从这一点上讲，编程与数学可谓关系密切。编程所需要的很多能力和数学是相通的，比如逻辑思维、抽象能力等。编程可以帮助程序员更好地理解数学，将复杂的数学可视化，也可以作为解决数学问题的工具，更能强化数学能力。本文作为一个系列的开篇之作，尝试站在程序员的角度，用程序员的方式去理解数学王国的经典概念。

如果有人问2的算术平方根是多少，相信所有的程序员都会张口说出1.414这个答案。不过，要是精度要求更高一点，大多数的程序员就只能依赖计算器了。比如，Python程序员会这样写：

>>> pow(2, 1/2)
1.4142135623730951

或者，像下面这样写，也没有问题：

>>> 2**(1/2)
1.4142135623730951

但是，如果想要通过计算的方式（不是指在草纸上手工开方），得到高精度的结果，作为程序员应该如何写代码呢？这里，以计算2的算术平方根为例，介绍科学巨人牛顿的逼近法。

牛顿逼近法是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法，虽然证明起来有点麻烦，但原理很简单，实用性很强。假设有方程：
$$f(x)=0$$

首先，选择一个接近函数$f(x)$零点的$x_0$，计算对应的$f(x_0)$和切线斜率$f^{{}^{\prime }}(x_0)$（这里$f^{{}^{\prime }}$表示函数$f$的
的导函数）；然后计算函数$f(x)$在点$(x_0,f(x_0))$的切线（斜率为$f^{{}^{\prime }}(x_0)$）和$x$轴的交点的$x$坐标，也就是求如下方程的解：
$$(x-x_0)f^{{}^{\prime }}(x_0)+f(x_0)=0$$
$$x=x_0-\frac{f(x_0)}{f^{{}^{\prime }}(x_0)}$$

将新求得的点的$x$坐标命名为$x_1$，通常$x_1$会比$x_0$更接近方程$f(x)=0$的解。因此可以利用$x_1$开始下一轮迭代。迭代公式可化简为如下所示：
$$x_n=x_{n-1}-\frac{f(x_{n-1})}{f^{{}^{\prime }}(x_{n-1})}$$

经过$n$次迭代，当$f(x_n)$的绝对值满足精度要求时，即可认为$x_n$就是方程$f(x)=0$的近似解。

理解了牛顿迭代法，就不难写出计算2的算术平方根的代码。

def f(x):
    # 定义函数f(x)
    
    return x**2 - 2

def df(x):
    # 定义函数f(x)的导函数
    
    return 2*x

def newton_raphson_method():
    # 牛顿-拉弗森方法
    
    tiny = 1e-15 # 当f(x_n)小于tiny时，迭代结束
    i, x = 0, 2 # 迭代计数器和初始值
    while True:
        i += 1 # 迭代计数器加1
        x = x-f(x)/df(x) # 计算下一个x
        if abs(f(x)) < tiny: # 满足迭代结束条件
            print('\\n经过%d次迭代，2的算术平方根为%.030f'%(i, x))
            break

if __name__ == '__main__':
    newton_raphson_method()

运行结果如下：

经过5次迭代，2的算术平方根为1.414213562373095145474621858739

请按任意键继续. . .

这个结果和使用pow函数计算出来的结果，小数部分的前30位完全一致。

>>> '%.030f'%pow(2,1/2)
'1.414213562373095145474621858739'

最后，再介绍一个计算函数导数的通用方法。对于类似$x^2-2$这样的函数，我们当然可以直接写出导函数$2x$。但如果函数足够复杂，不能直接写出导函数，则可以象下面这样使用微分法近似计算函数在$x=x_i$处的导数。

def df(xi):
    # 返回函数f(x)在xi处的导数
    
    delta = 1e-15
    return (f(xi+delta)-f(xi-delta))/(2*delta)

使用这个近似的求导函数，迭代得到的结果就会有一点误差，小数部分的前15位和上面两种方法一致。

经过15次迭代，2的算术平方根为1.414213562373095367519226783770

请按任意键继续. . .