Skip to content

Latest commit

 

History

History
1161 lines (741 loc) · 39.2 KB

File metadata and controls

1161 lines (741 loc) · 39.2 KB

十一、Python 的更多微积分

概述

在本章中,您将学习如何计算给定方程式的曲线长度。将向您介绍三维中的偏导数,以及如何使用它们来计算曲面的面积。跟随中世纪数学家的脚步,您将使用一个无穷级数来计算常数,如 pi,并确定级数的收敛间隔。就像现代数学家和机器学习工程师一样,您将学习如何使用偏导数找到曲面上的最小点。到本章结束时,您将能够使用微积分来解决各种数学问题。

导言

在上一章中,我们学习了如何计算导数和积分。现在,我们将使用这些工具来找到曲线和螺旋的长度,并将这种推理扩展到三维,以找到复杂曲面的面积。我们还将介绍微积分中使用的一个常用工具,无穷级数,它用于计算重要常数和近似复杂函数。最后,我们将研究机器学习中的一个重要思想:找到曲线上的最小点。当你使用神经网络时,你会创建一种“误差函数”,努力寻找表面上误差最小的点。我们将创建我们自己的梯度下降函数,继续向下移动,直到到达曲面底部。

曲线的长度

导数和积分的主要用途是求曲线的长度。这有一个公式:

Figure 11.1: Formula to calculate the length of a curve

图 11.1:计算曲线长度的公式

前面的公式包含一个积分和一个导数。为了求曲线的长度,我们需要导数和积分函数。如果尚未将它们复制并粘贴到代码中:

from math import sqrt

def derivative(f,x):
    """Returns the value of the derivative of     the function at a given x-value."""
    delta_x = 1/1000000
    return (f(x+delta_x) - f(x))/delta_x

def trap_integral(f,a,b,num):
    """Returns the sum of num trapezoids     under f between a and b"""
    width = (b-a)/num
    area = 0.5*width*(f(a) + f(b) + 2*sum([f(a+width*n) \                                        for n in range(num)]))
    return area

下面是公式的 Python 版本:

def curve_length(f,a,b,num):
    def g(x):
        return sqrt(1+(derivative(f,x)**2))
    return trap_integral(g,a,b,num)

注意,我们只是将数学符号转换为 Python 代码。我们在f函数中定义了g函数。g函数是公式中平方根下的一切。然后,我们使用我们的trap_integral函数找到ab之间的g函数的累积值。

让我们用一条我们知道长度的曲线来检查,比如直线y=2x。我们可以使用毕达哥拉斯定理计算曲线的直线在*x=(0,0)x=(2,4)*之间的距离√5 或 4.47 单位:

def f(x):
    return 2*x
print(curve_length(f,0,2,1000))

前面的代码打印出 4.47。。。作为输出。

但是,当我们试图检查一条我们知道长度的实际曲线时,例如半圆,我们遇到了一个问题。我们知道下式曲线的长度,因为它是半径为 1 的圆周长的一半。所以,我们应该得到π或 3.1415。。。作为输出:

Figure 11.2: Formula to calculate the length of a semicircle

图 11.2:计算半圆长度的公式

让我们把f(x)换成前面半圆形的方程式:

def f(x):
    return sqrt(1-x**2)
print(curve_length(f,-1,1,100))

当您执行前面的代码时,您会得到一个错误。错误消息的最后一行(我读到的第一行)说:

ValueError: math domain error

这是因为-1 和 1 处半圆的导数是无穷大的。这些点的切线是垂直的,如下图所示:

Figure 11.3: Vertical tangent lines, with an infinite slope

图 11.3:具有无限斜率的垂直切线

因此,这种方法已经遇到了问题。让我们看看它是否能找到正则多项式的长度,如下图所示:

Figure 11.4: A complicated polynomial

图 11.4:复杂多项式

这是一个 5 次多项式,表示x的最高指数为5。曲线方程如下所示:

Figure 11.5: Equation of the curve

图 11.5:曲线方程

尽管看起来很复杂,但曲线上没有一个导数是无限的,正如图 11.3所示。这意味着我们可以在上面使用曲线长度函数。

以下是多项式的代码:

def f(x):
    return 0.7*x**5 + 1.6*x**4-2.05*x**3 -3*x**2+2.95*x+2.9
print(curve_length(f,-2,1,1000))

曲线的长度如下所示:

9.628984854276812

我们可以使用 Wolfram Alpha 来解决这个问题,方法是输入曲线 y 的长度=。。。从-2 到 1并检查它是否是一个良好的近似值。但是使用 Python,有一种更直接的计算曲线长度的方法,它不会遇到我们在使用导数时遇到的问题。事实上,它甚至不使用导数或积分。您可以使用勾股定理简单地求出曲线的一小段长度,然后将所有这些小段相加,如下图所示。我们知道宽度,我们对小直角三角形的斜边感兴趣。我们可以计算高度,即x处的函数与x处的函数之差,加上宽度:

Figure 11.6: Finding the length of a tiny part of a curve

图 11.6:找到曲线一小部分的长度

上图所示的直角三角形斜边如下:

图 11.7:计算直角三角形斜边的公式

我们所要做的就是通过从ab的间隔,计算所有这些长度。以下是如何在 Python 中实现这一点:

def f(x):
    return 0.7*x**5 + 1.6*x**4-2.05*x**3 -3*x**2+2.95*x+2.9
def curve_length2(f,a,b,num=1000):
    """Returns the length of f between\
    a and b using num slices"""
    output = 0
    width = (b-a)/num
    for i in range(num):
        output += sqrt((f(a+(i+1)*width)-f(a+i*width))**2 + width**2)
    return output

这应该会提醒您积分程序:创建一个连续的和,然后在曲线的每个切片上循环,同时添加面积(在本例中为弧长)。最后,返回运行总和的最终值。

下面是我们感兴趣的区间曲线长度:

print(curve_length2(f,-2,1))

这给了我们曲线的长度为9.614118659973549。这比以前的版本更接近,并且没有太多的麻烦。现在,轮到你在下面的练习中做同样的事情了。

练习 11.01:寻找曲线的长度

在本练习中,将为您提供以下曲线方程式。使用该方程,确定两个给定的x值之间的曲线长度:

Figure 11.8: Equation of the curve

图 11.8:曲线方程

这些值将从x=-1x=1

执行以下步骤:

  1. First, we need to create a circle function with the preceding equation:

    def circle(x):
        return sqrt(1-x**2)

    笔记

    又是半圆了。这一次,我们的curve_length2函数将不会有任何问题,将微小的弧片相加。

  2. 现在,我们将在该曲线上运行curve_length2函数(我们已经对其进行了编码),将所有微小的线段相加,正如我们之前所做的:

    def curve_length2(f,a,b,num=1000):
        """Returns the length of f between\
           a and b using num slices"""
        output = 0
        width = (b-a)/num
        for i in range(num):
            output += sqrt((f(a+(i+1)*width)-f(a+i*width))**2 \
                            + width**2)
        return output
  3. Now, we print the output of the function, measuring from x = -1 to x = 1:

    print(curve_length2 (circle,-1,1))

    结果如下:

    3.1415663562164773

这次没有错误消息。我们得到半径为 1 的圆的半圆周长的一个很好的近似值,我们知道它是π。

笔记

要访问此特定部分的源代码,请参考https://packt.live/3gkI5Qi

您也可以在在线运行此示例 https://packt.live/3eVpSbz

练习 11.02:寻找正弦波的长度

正弦波是数学和科学中一个非常重要和有用的函数。它在 0 和 2π之间进行一个循环,如下图所示:

Figure 11.9: One cycle of the sine wave

图 11.9:正弦波的一个周期

很容易测量它的波长(2π)和振幅(它向上和向下走多远,即 1 个单位),但实际曲线有多长?在本练习中,我们将找到从 0 到 2π的正弦波长度。

执行以下步骤:

  1. 我们将再次使用我们的curve_length2函数,但现在我们必须从math模块

    from math import sin, pi

    导入我们的sinpi函数

  2. We've already written our curve_length2 function, which will add up all the segments of the curve. We just need to tell it the function to use, and the beginning and ending x values:

    print(curve_length2(sin,0,2*pi))

    结果如下:

    7.640391636335927

如您所见,使用curve_length2函数,计算正弦波的长度变得非常容易。

笔记

要访问此特定部分的源代码,请参考https://packt.live/3dUy3nk

您也可以在在线运行此示例 https://packt.live/2VFy2xd

螺旋的长度

那么用极坐标表示的螺旋呢,其中r是到原点的距离,是与x轴形成的θ(θ)角的函数?我们不能使用我们的xy函数来测量下图所示的螺旋线:

Figure 11.10: An Archimedean spiral

图 11.10:阿基米德螺线

我们在上图中看到的是一个螺旋,从(5,0)开始,旋转 7.5 圈,以(11,π)结束。该曲线的公式为r(θ)=5+0.12892θ。旋转的弧度数是 2π的 7.5 倍,即 15π。我们将使用与上一节相同的想法:我们将找到从*r(θ)r(θ+步长)*的直线长度,用于中心角的某个小步,如下图所示:

Figure 11.11: Approximating the length of a tiny part of the curve

图 11.11:曲线一小部分的近似长度

上图所示三角形中心角的另一侧与积分问题中的切片或之前的曲线长度程序中三角形的斜边相似。这一次,它不是直角三角形,所以我们不能使用斜边。但对于这个问题,有一个公式叫做余弦定律。在三角形 ABC 中,对边角 C 的长度如下:

Figure 11.12: Law of cosines

图 11.12:余弦定律

我们需要做的就是把它放到一个函数中,像这样:

def opposite(a,b,C):
    """Returns the side opposite the given angle in
       a triangle using the Law of Cosines
       Enter side, side, angle"""
    c = sqrt(a**2 + b**2 - 2*a*b*cos(C))
    return c

然后,我们只需要写一个函数,从起始角度开始,沿着曲线走一小步,测量每个小角度对面的边,直到到达终止角度:

from math import sqrt,cos,pi
def spiral(r,a,b,step=0.0001):
    """Returns length of spiral r from
       a to b taking given step size"""
    length = 0
    theta = a
    while theta < b:
        length += opposite(r(theta),r(theta+step),step)
        theta += step
    return length

我们的职能如下:

def r(theta):
    return 5 + 0.12892*theta

我们所要做的就是在这个螺旋上执行螺旋函数,从 0 到 15π:

spiral(r,0,2*pi*7.5)

结果如下:

378.8146271783955

如您所见,螺旋的长度为378.8146271783955。在下一个练习中,我们将了解如何找到极轴缓和曲线的长度。

练习 11.03:找到极性螺旋曲线的长度

在本练习中,您将找到极轴缓和曲线的长度,该曲线从(3,0)开始,绕中心旋转 12 圈,到(16,0)结束。

执行以下步骤以查找所需的长度:

  1. We don't know the formula for this spiral, but we do know that the radius increases 13 units (from 3 to 16) in 12 revolutions. This means that for every increase of 2π in the angle, θ, the radius increases 13/12 units. So, we divide 13/12 by 2π. The increase in radius can be expressed as follows:

    Figure 11.13: Formula to calculate the increase in radius

    图 11.13:计算半径增加的公式

  2. 我们可以用 Python 这样表示:

    def r(theta):
        return 3 + 0.1724*theta
  3. We can check to make sure r(0) = 3 and r(24π) = 16 this way:

    print(r(0),r(24*pi))

    结果如下:

    3.0 15.998653763493127
  4. Now, we simply put that in our spiral function:

    spiral(r,0,2*pi*12)

    结果如下:

    716.3778471288748

在本练习中,我们很容易找到这条螺旋曲线的长度,即716.3778471288748,只需知道曲线的起点和终点值以及围绕中心旋转的圈数即可。

笔记

要访问此特定部分的源代码,请参考https://packt.live/2YT70EH

您也可以在在线运行此示例 https://packt.live/2YV4wFT

练习 11.04:查找卷中绝缘层的长度

已要求您找到下图所示的卷中剩余绝缘层的(近似)长度:

Figure 11.14: Measuring rolled up materials using calculus

图 11.14:使用微积分测量卷起材料

测量辊,发现中心是一个直径为 4 英寸的空圆(因此r(0)=2。辊的外径为 26 英寸。你从中心向外数层,估计螺旋需要 23.5 圈,因此r(2π23.5)=26/2=13*。

执行以下步骤以计算长度:

  1. Calculate the equation using the preceding data:

    Figure 11.15: Formula to calculate radius

    图 11.15:计算半径的公式

  2. Here's what the graph of the spiral looks like:

    Figure 11.16: A graph of the roll of insulation

    图 11.16:绝缘辊图

  3. 现在,不难将我们的r代码更改为这个螺旋:

    def r(theta):
        return 2 + 0.0745*theta
  4. Now, we can run our spiral function on this function from 0 to 2π23.5:

    spiral(r,0,2*pi*23.5)

    结果如下:

    1107.502879450013

1107 英寸的绝缘层刚刚超过 92 英尺。

笔记

要访问此特定部分的源代码,请参考https://packt.live/2VE9YKZ

您也可以在在线运行此示例 https://packt.live/31D43tG

练习 11.05:寻找阿基米德螺旋的长度

在这个练习中,你已经得到了阿基米德螺线的方程。从θ=0θ=2π的螺旋线长度:

Figure 11.17: Equation of an Archimedean spiral

图 11.17:阿基米德螺线方程

笔记

这适用于对数螺线和阿基米德螺线。

执行以下步骤以查找长度:

  1. 我们简单地用指数函数重新定义了r

    from math import e
    def r(theta):
        return 2*e**(0.315*theta)
  2. Then. we run the spiral function from 0 to :

    spiral(r,0,2*pi)

    结果如下:

    41.518256747758976

该螺旋的长度为41.518256747758976

笔记

要访问此特定部分的源代码,请参考https://packt.live/2VEtjfo

您也可以在在线运行此示例 https://packt.live/2VHasQN

表面的面积

让我们学习如何从二维到三维,并计算三维曲面的面积。在第 10 章**基础微积分与 Python中,我们学习了如何计算旋转曲面的面积,但这是一个曲面,其中三维zxy值的函数。

公式

通过曲面上的二重积分给出了解析求解该问题的传统代数方法:

Figure 11.18: Formula to calculate area of a surface

图 11.18:计算表面面积的公式

这里,z=f(x,y)(x,y,f(x,y))。这些卷曲的 d 是三角洲,这意味着我们将处理偏导数。偏导数是对一个变量的导数,即使函数依赖于多个变量。这是一个函数,它返回函数f在特定点(v,w)相对于变量u的偏导数。根据我们感兴趣的变量,xy,函数将朝该方向迈出一小步并计算导数,正如我们已经做过的:

def partial_d(f,u,v,w,num=10000):
    """returns the partial derivative of f
    with respect to u at (v,w)"""
    delta_u = 1/num
    try:
        if u == 'x':
            return (f(v+delta_u,w) - f(v,w))/delta_u
        else:
            return (f(v,w+delta_u) - f(v,w))/delta_u
    except ValueError:
        pass

如果抛出一个ValueError,代码中有一个try...except块。如果坡度太大,就会发生这种情况,如在垂直线上。如果发生这种情况,它将忽略它并继续。

现在,我们需要一个 3D 向量和一个cross函数,用于面积公式中的叉积。叉积给出垂直于两个给定向量的向量长度,以及由给定向量形成的平行四边形的面积:

Figure 11.19: The cross product of two vectors

图 11.19:两个向量的叉积

如果知道向量之间的角度,可以使用该角度来查找叉积:

Figure 11.20: Formula to calculate cross product of two vectors

图 11.20:计算两个向量的叉积的公式

如果没有,就像我们的例子一样,可以使用 3D 向量来表示向量在各个方向上的位移,xyz。例如,假设我们有两个向量,u=2i+3j+4kv=5i+6j+7k。它们由它们在三个维度中的位移来定义。i部分为x方向的位移;j部分为y方向的位移;k部分为z方向的位移。好消息是会有一些零来简化事情。要交叉两个向量,我们可以将它们放入一个矩阵中,并将它们作为下列矩阵的行列式进行运算:

Figure 11.21: Calculating cross product of two vectors using matrix

图 11.21:使用矩阵计算两个向量的叉积

我们将编写一个函数来对两个 3D 向量执行该操作。我们需要输入的只是ijk的系数。因此,如果u=ai+bj+ckv=di+ej+fk,我们将得到以下结果:

Figure 11.22: Performing mathematical operations on 3D vectors

图 11.22:对 3D 向量执行数学运算

让我们使用列表来表示向量,因此系数为u=[a,b,c]u[0]=au[1]=bu[2]=c

def cross(u,v):
    """Returns the cross product of 2 3D vectors
    [[i,j,k],\
    [1,0,dz/dx],\
    [0,1,dz,dy]]
    cross([1,-1,2],[2,3,-5])
    >>> [-1, -9, 5]
    """
    return [u[1]*v[2]-v[1]*u[2],\
            -u[0]*v[2]+v[0]*u[2],\
            u[0]*v[1]-v[0]*u[1]]

我们编写了一个长 docstring,以明确函数的用途、如何输入值以及我们将得到什么作为输出。让我们检查一下,以确保得到正确的输出:

print(cross([2,3,4],[5,6,7]))

结果如下:

[-3, 6, -3]

这很有效。现在,我们需要写一个函数来求一个 3D 向量的大小,因为这将给出平行四边形的面积。这只是毕达哥拉斯定理在三维空间的延伸。因此,如果u=ai+bj+ck的向量u的大小为a

def mag(vec):
    """Returns the magnitude of a 3D vector"""
    return sqrt(vec[0]**2+vec[1]**2+vec[2]**2)

这是半圆的样子,它的表面近似于平行四边形。更多平行四边形意味着更精确的近似值:

Figure 11.23: Using more parallelograms

图 11.23:使用更多平行四边形

我们的面积函数将循环通过网格中的所有xy点,计算每个点的偏导数,并使用叉积计算该点平行四边形的面积:

from math import sqrt
def sphere(x,y):
    """Sphere of radius 1"""
    return sqrt(1-x**2-y**2) 
def area(f,ax,bx,ay,by,num=1000):
    """Returns area of parallelogram formed by
    vectors with given partial derives"""
    running_sum = 0
    dx = (bx-ax)/num
    dy = (by-ay)/num
    for i in range(num):
        for j in range(num):
            x = ax+i*dx
            y = ay+j*dy
            dz_dx=partial_d(f,'x',x,y)
            dz_dy=partial_d(f,'y',x,y)
            try:
                running_sum += mag(cross([1,0,dz_dx],[0,1,dz_dy]))*dx*dy
            except:
                pass
    return running_sum

首先,我们将区域的运行和设置为 0。然后,我们计算dxdyxy中的微小变化,将曲面分成相等的切片。当与球体相切的直线的斜率垂直时,try...except块简单地忽略(pass)偏导数为无穷大时将产生的误差,如图 11.3所示。如果没有错误,它会加上在该点由偏导数形成的平行四边形的面积。现在,我们在半球上运行面积函数,在每个方向上使用 1000 个点,得到一个非常精确的近似值。我们知道半径为 1 的球体的一半表面积是 2π,或 6.28:

print("Area of hemisphere:",area(sphere,-1,1,-1,1))

结果如下:

Area of hemisphere: 6.210356913122

现在,让我们快速执行基于此概念的练习。

练习 11.06:寻找 3D 曲面的面积-第 1 部分

现在,让我们找到一个复杂曲面的面积,这很难用代数方法找到。考虑下列方程的表面:

Figure 11.24: Equation of the surface

图 11.24:表面方程

该曲面如下图所示:

Figure 11.25: A complicated 3D surface

图 11.25:复杂的 3D 曲面

执行以下步骤以查找该区域:

  1. 让我们把这个函数放到我们的区域程序中,看看我们得到了什么:

    from math import sin, cos, sqrt
    def surface(x,y):
        return 10*sin(sqrt(x**2+y**2))
    print("Area of wave surface:",area(surface,-5,5,-5,5))
  2. 运行程序查看输出:

    Area of wave surface: 608\. 2832236305994

查看前面的代码,我们可以清楚地看到,使用 Python 只需几行代码就可以很容易地找到复杂曲面的区域。

笔记

要访问此特定部分的源代码,请参考https://packt.live/3gwd6kr

您也可以在在线运行此示例 https://packt.live/2ZpgwOQ

练习 11.07:寻找 3D 曲面的面积-第 2 部分

找到曲面的面积a

下面是曲面的外观:

Figure 11.26: Another 3D surface

图 11.26:另一个 3D 曲面

执行以下步骤以查找该区域:

  1. 定义曲面函数以返回表达式:

    def surface(x,y):
        return 3*cos(x)+2*cos(x)*cos(y) 
  2. Run the surface function to get the value:

    print("Area of surface:",area(surface,0,6.28,0,6.28))

    结果如下:

    Area of surface: 99.80676808568984

该 3D 曲面的面积为99.80676808568984

笔记

要访问此特定部分的源代码,请参考https://packt.live/2VCaObq

您也可以在在线运行此示例 https://packt.live/2NPXvQo

练习 11.08:寻找表面区域–第 3 部分

找到曲面的面积b

下面是曲面的外观:

Figure 11.27: The surface of

图 11.27a的表面

执行以下步骤以查找该区域:

  1. 定义曲面函数以返回新表达式:

    def surface(x,y):
        return sqrt(1+sin(x)*cos(y))
  2. Run the surface function:

    print("Area of surface:",area(surface,0,6.28,0,6.28))

    结果如下:

    Area of surface: 42.80527549685105

该表面的面积为42.80527549685105

笔记

要访问此特定部分的源代码,请参考https://packt.live/3gwdLlV

您也可以在在线运行此示例 https://packt.live/3dUNWdt

无穷级数

数学家经常遇到太复杂的函数,无法求解或以其他方式处理,而近似一直是数学中的一个重要组成部分。对于试图用代数方法求导数和积分的数学家来说,许多表达式都没有很好的精确解、导数、积分等等。一般来说,科学家在现实生活中遇到的微分方程都没有代数解,所以他们必须使用其他方法。稍后将详细介绍微分方程,但有一个重要的近似家族使用简单函数来近似函数。

多项式函数

很容易求解、微分和积分多项式方程,例如y=x2,甚至以下方程:

Figure 11.28: A polynomial equation

图 11.28:多项式方程

这些项都是一个接一个地加(或减)的,并且没有三角函数、对数函数或指数函数使事情变得困难。以下是用简单多项式逼近函数的公式:

Figure 11.29: The Taylor series

图 11.29:泰勒级数

这个公式被命名为泰勒级数:任何函数(可导函数)都可以在特定点上以一定精度近似,只需使用具有特定页面的多项式即可。

系列

数学家有一个符号来表示将一组遵循某种模式的数字相加:

Figure 11.30: Formula for calculating series

图 11.30:计算系列的公式

看起来像E的大符号实际上是希腊字母 sigma,或S,它表示数字的总和。西格玛下方的等式是变量开始的位置(在本例中为 1),而西格玛上方的等式是i(在本例中为 10)的最后一个整数值。西格玛右边是一个表达式,表示如何处理变量。在本例中,我们只是添加变量i,从 1 到 10。这几乎就是用 Python 编写列表理解的方式。以下是方法:

s = sum([i for i in range(1,11)])

列表理解中的第一个术语是您在本例中的西格玛级数表达式中看到的,i。例如,整数的平方和到n的级数如下:

Figure 11.31: Series for sum of squares of integers from 1 to n

图 11.31:从 1 到 n 的整数平方和序列

在 Python 中,我们将这样编写:

s = sum([i**2 for i in range(1,n+1)])

一个古老但有用的数列是反正切数列。它计算具有给定切线的角度(以弧度为单位);例如:

Figure 11.32: Equation of an arctangent series

图 11.32:反正切级数的方程

根据上述方程式,arctan 的方程式如下:

Figure 11.33: Equation of an arctan

图 11.33:arctan 方程

该系列由以下模式计算:

Figure 11.34: Equation for series of arctan

图 11.34:arctan 系列方程

下面是 sigma 表达式:

Figure 11.35: A sigma expression

图 11.35:σ表达式

通过插入x的切线,我们可以计算出角度的近似值:

Figure 11.36: Substituting the value of x in the equation

图 11.36:替换方程式中的 x 值

这对于几个世纪前的数学家来说是一个相当大的计算量,但这里是 Python 的等价物。注意列表理解的第一部分与前面的 sigma 表达式有多接近:

def arctan(x,n):
    """Returns the arctangent of x using a series of n terms."""
    return sum([((-1)**(i-1)*(x**(2*i-1)))/(2*i-1) \
                              for i in range(1,n+1)])
print(arctan(1/1.732,10))

因此,在 10 个条款之后,我们得到以下结果:

0.523611120446175

这与a非常接近。

收敛

数学家希望简化 arctan 级数,以便使用以下事实轻松计算π

Figure 11.37: Trigonometric function of a tangent

图 11.37:切线的三角函数

根据上述方程式,arctan 的方程式如下:

Figure 11.38: Formula to calculate arctan

图 11.38:计算 arctan 的公式

他们认为在 arctan 系列中用1替换x将使计算π成为公园里的散步。以下是前几个术语:

Figure 11.39: Substituting x = 1 in the arctan series

图 11.39:在 arctan 系列中替换 x=1

该表达式给出了pi的近似值:

Figure 11.40: Equation to find the approximate value of pi

图 11.40:求 pi 近似值的方程式

我们只需在 sigma 右边写零件代码,在n范围内添加代码,并将其汇总:

for n in range(1,10):
    print(4*sum([((-1)**(i-1))/(2*i-1) for i in range(1,n+1)]))

我们可以在输出中显示接近pi的进度:

4.0
2.666666666666667
3.466666666666667
2.8952380952380956
3.3396825396825403
2.9760461760461765
3.2837384837384844
3.017071817071818
3.2523659347188767

这与π不是很接近。跳过更大数量的术语如何?让我们更改循环的代码:

for n in [100,1000,1000000]:

这是输出:

3.1315929035585537
3.140592653839794
3.1415916535897743

100 万个术语之后,我们只得到了正确的小数点后五位。该系列收敛到(即非常接近或趋向)π/4 的速度对于任何实际用途来说都太慢。几个世纪以来,数学家们一直在寻找更好的数列来近似π。

练习 11.09:计算π的 10 个正确数字

1706 年,英国数学家和天文学家约翰·梅切尔(John Machine)用他的改进数列计算了π的 100 位小数。以下是系列:

Figure 11.41: An arctan function

图 11.41:arctan 函数

使用前面的 arctan 函数计算 10 位正确的π。请按照以下步骤执行此操作:

  1. 只需调用我们的 arctan 函数。10 个条款应足够:

    print(4*(4*arctan(1/5,10)-arctan(1/239,10)))
  2. 运行前面的代码查看输出:

    3.1415926535897922

我们用 10 个项得到了一个很好的近似值。它给出的正确数字甚至超过 10 个。

笔记

要访问此特定部分的源代码,请参考https://packt.live/3dPjVvD

您也可以在在线运行此示例 https://packt.live/3dVlTKR

练习 11.10:使用欧拉表达式计算π的值

伟大的德国数学家欧拉提出了以下表达式:

Figure 11.42: Euler's expression

图 11.42:欧拉表达式

使用此表达式近似于π。它是否比调整后的 arctan 公式收敛更快?

执行以下步骤:

  1. 以下是使用欧拉级数逼近π的代码:

    from math import sqrt
    for n in [100,1000,1000000]:
        print(sqrt(6*sum([1/(i**2) for i in range(1,n+1)])))
  2. 它会聚得更快吗?运行前面的代码查看输出:

    3.1320765318091053
    3.1406380562059946
    3.1415916986605086

它似乎并没有更快地收敛。在 100 万个术语之后,您仍然只有五个正确的小数位。

笔记

要访问此特定部分的源代码,请参考https://packt.live/2NRnnLD

您也可以在在线运行此示例 https://packt.live/38lHXgm

20 世纪的公式

下面是杰出的印度数学家 Ramanujan 的π近似公式:

Figure 11.43: Ramanujan's formula to approximate π

图 11.43:Ramanujan 近似π的公式

以下是如何用 Python 编写代码:

from math import sqrt, factorial
one_over_pi = 2*sqrt(2)/9801*sum([(factorial(4*k)*(1103+26390*k))/ \
              (((factorial(k))**4)*(396**(4*k))) for k in range(10)])
print(1/one_over_pi)

10 项后的输出如下:

3.141592653589793

那很准确!!

收敛区间

级数收敛(趋向于某个值)的值范围称为收敛区间。使用 Python,找到这个区间相当简单:在序列中运行一些数字,如果它们变得无限大,它们就不在区间内。如果它们产生一个数字,它们就在区间内。例如,让我们看一个非常常见的教科书问题,并使用 Python 解决它。

练习 11.11:确定收敛间隔-第 1 部分

确定以下幂级数的收敛间隔:

Figure 11.44: A power series

图 11.44:幂级数

执行以下步骤:

  1. Enter the sum into Python:

      def mystery_sum(x):
        return sum([(((-1)**n)*n)/(4*n)*(x+3)**n for n in \
                          range(1,1000000)])

    因为我们不能使用“无穷大”这个数字,所以我们可以找到 n=1 到 100 万之间的所有项的总和。

  2. 运行-10 到 10 之间的所有整数,查看是否有收敛到某个数字的整数:

    for x in range(-10,11):
        print(x,mystery_sum(x))
  3. When you run this, you'll get an OverflowError:

    OverflowError: int too large to convert to float

    所有这一切意味着一些数字变得越来越大,这是我们所期望的。我们需要添加一个条件,这样如果我们得到那个错误,它将简单地返回Infinity。这是通过try...except块完成的。

  4. 让我们告诉 Python 尝试一行代码。如果它抛出一个特定错误(在本例中为OverflowError),不要停止程序,只需执行以下操作:

    def mystery_sum(x):
        try:
            return sum([(((-1)**n)*n)/(4*n)*(x+3)**n \
                                      for n in range(1,1000000)])
        except OverflowError:
            return "Infinity"
  5. 现在,输出给我们一些无穷大和一些实际值:

    -10 Infinity
    -9 Infinity
    -8 Infinity
    -7 Infinity
    -6 Infinity
    -5 Infinity
    -4 249999.75
    -3 0.0
    -2 -0.25
    -1 Infinity
    0 Infinity
    1 Infinity
    ...

看起来我们的收敛区间是-5< x < -1. This means we can use the series to get useful values if x在区间内。否则,我们不能使用它。

笔记

要访问此特定部分的源代码,请参考https://packt.live/38k30A2

您也可以在在线运行此示例 https://packt.live/31AtmMU

练习 11.12:确定收敛间隔-第 2 部分

确定以下幂级数的收敛间隔:

Figure 11.45: A power series

图 11.45:幂级数

执行以下步骤:

  1. 用 Python 定义总和:

    def mystery_sum(x):
        try:
            return sum([n*x**n/(5**(2*n)) for n in range(1,10000)])
        except OverflowError:
            return "Infinity"
    
    for x in range(-30,30):
        print(x,mystery_sum(x))
  2. 以下是一些输出:

    -30 Infinity 
    -29 Infinity 
    -28 Infinity 
    -27 Infinity 
    -26 -1.0561866634327267e+174 
    -25 -5000.0 
    -24 -0.24989587671803576 
    -23 -0.24956597222222246 
    -22 -0.24898143956541371 
    -21 -0.24810964083175827 
    -20 -0.24691358024691298
    ...
    18 9.18367346938776 
    19 13.19444444444444 
    20 19.999999999999993 
    21 32.812499999999964 
    22 61.11111111111108 
    23 143.74999999999983 
    24 599.9999999999994 
    25 49995000.0 
    26 5.3728208568640556e+175 
    27 Infinity 
    28 Infinity
    29 Infinity

在-25 和 25 之间的x的所有输出都很小(在 0 和 600 之间),无论我们使用了多少项,因此我们将收敛间隔称为*-25<x<25*。

笔记

要访问此特定部分的源代码,请参考https://packt.live/38pwuwC

您也可以在在线运行此示例 https://packt.live/2YS46jl

练习 11.13:寻找常数

在本练习中,我们将用 Python 表示一个无穷级数并求和。我们将使用一个著名的常数,它被定义为系列的总和:

Figure 11.46: Sum of the series

图 11.46:系列总和

这个著名常数的值是多少?让我们按照以下步骤确定该值:

  1. 导入阶乘模块,将前面的公式转换成 Python,如下:

    from math import factorial
    print(sum([1/factorial(n) for n in range(10000)]))
  2. 运行前面的代码查看输出:

    2.7182818284590455

著名的常数是e,自然对数的基数。

笔记

要访问此特定部分的源代码,请参考https://packt.live/2AoyubH

您也可以在在线运行此示例 https://packt.live/2BZ4aVw

活动 11.01:找到表面的最小值

机器学习的一个主要任务是最小化函数。当你训练神经网络时,你要改变矩阵或张量中的值,看看哪一个能更好地逼近你的测试数据。在网络中的每个值上,您都可以看到它对错误的贡献程度。这听起来像是曲面上不同点的偏导数,不是吗?

梯度下降过程就是一个例子。让我们考虑一下,我们要找到函数的最小值的位置。曲面上的每一点都有偏导数,我们可以利用偏导数稍微向一个较低的值移动。我们将从随机的某个地方开始,计算该点的偏导数,然后沿着降低z值的方向移动,即上下值。所以,如果zx(我们称之为dz_dx的偏导数为负,这意味着z随着x的增加而减少,我们希望向正x方向移动。如果dz_dx为正,这意味着z随着x的增加而增加,所以我们要朝相反的方向移动,所以我们将朝负x方向移动。我们将对y方向执行同样的操作。这将如下所示:

Figure 11.47: The path of descent to a minimum value

图 11.47:下降至最小值的路径

此活动的第一部分是创建一个函数,用于查找曲面的最小点。可通过以下步骤编写此函数:

  1. 编写一个函数,在曲面上创建一个随机(x,y)位置。您可以调用random模块的uniform函数来生成这些值。
  2. 计算zxy的偏导数。
  3. xy更改为偏导数的负数,如果偏导数较大,则乘以微小的步数
  4. 在此新位置计算偏导数,并保持循环,直到偏导数都非常小(小于 0.0001)或该位置离开曲面。
  5. 在一组随机位置上运行该函数,将最小的z值保存到分钟列表中。
  6. 最后,打印最小的分钟列表。

编写函数后,请在已知函数值的函数上进行测试,以验证其是否按预期工作。然后,您可以在不知道最小点的函数上运行它,以确定此未知位置。具体步骤如下:

  1. 在地面6上测试您的功能。您的函数应该在点(0,0)处发现最小值为 0。
  2. 一旦你对自己的功能有信心,用它来确定*-1<x<5-1<y<5*的最小值。

你会发现,根据你的起点,你的函数会收敛到不同的最小值点——局部最小值和全局最小值。

笔记

可在第 696 页找到此活动的解决方案。

总结

在上一章中,我们学习了导数和积分的威力,因此在本章中,我们利用这些工具来解决一些相当困难的问题,例如螺旋的长度和 3D 曲面的面积。我们甚至通过引入偏导数将导数和积分推广到三维。在微积分类中,为了使用这些工具,我们将使用大量代数,但通过使用 Python,我们对情况进行了建模并测试了我们的函数。我们创建了包含不断变化的值的变量,并在必要时循环计算数百万次。对于前几个世纪的数学家来说,这似乎是一种神奇的灯。

在下一章中,我们将处理更多的变化率和数量,并通过使用 Python 避免大量代数运算。我们将了解在一种不断变化的混合物中有多少盐,捕食者何时何地捕捉猎物,以及我们需要投资多长时间才能赚到 100 万美元。

FAB62

RUC47