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

曾经终败给现在 2023-01-18 12:20 39阅读 0赞

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

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

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

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

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

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

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

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

将新求得的点的 x x x坐标命名为 x 1 x\_1 x1，通常 x 1 x\_1 x1会比 x 0 x\_0 x0更接近方程 f ( x ) = 0 f(x)=0 f(x)=0的解。因此可以利用 x 1 x\_1 x1开始下一轮迭代。迭代公式可化简为如下所示：  
 x n = x n − 1 − f ( x n − 1 ) f ′ ( x n − 1 ) x\_n = x\_\{n-1\} - \\frac\{f(x\_\{n-1\})\}\{f^\{'\}(x\_\{n-1\})\} xn=xn−1−f′(xn−1)f(xn−1)

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center]  
理解了牛顿迭代法，就不难写出计算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 x^2-2 x2−2这样的函数，我们当然可以直接写出导函数 2 x 2x 2x。但如果函数足够复杂，不能直接写出导函数，则可以象下面这样使用微分法近似计算函数在 x = x i x=x\_i x=xi处的导数。

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

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

经过15次迭代，2的算术平方根为1.414213562373095367519226783770
    
    请按任意键继续. . .

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center]: /images/20221021/06411d06a69a4452ae38e818163fbb7a.png