Python-统计和微积分研讨会（六）

Python 统计和微积分研讨会（六）

原文：zh.annas-archive.org/md5/6cbaed7d834977b8ea96cc7aa6d8a083

译者：飞龙

协议：CC BY-NC-SA 4.0

第十一章：11.使用 Python 进行更多微积分

概述

在本章中，您将学习如何根据其方程计算曲线的长度。您将学习如何在三维空间中使用偏导数来计算表面积。跟随中世纪的数学家的脚步，您将使用无限级数来计算诸如π之类的常数，并确定级数的收敛区间。像现代数学家和机器学习工程师一样，您将学习如何使用偏导数找到表面上的最小点。在本章结束时，您将能够使用微积分解决各种数学问题。

介绍

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

曲线的长度

导数和积分的一个主要用途是找到曲线的长度。有一个公式：

图 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函数来找到a和b之间g函数的累积值。

让我们检查一下我们知道长度的曲线，比如线y = 2x。我们可以使用勾股定理计算曲线的长度，x = (0,0)和x = (2,4)之间的距离，结果为 2√5 或 4.47 个单位：

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

前面的代码输出了 4.47...作为输出。

但是当我们尝试检查我们知道长度的实际曲线，比如半圆时，我们遇到了问题。我们知道以下方程的长度，因为它是半径为 1 的圆的一半周长。所以，我们应该得到π或 3.1415...作为输出：

图 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 处的导数是无穷大。这些点处的切线是垂直的，如下图所示：

图 11.3：垂直切线，斜率无穷大

所以，这种方法已经遇到了问题。让我们看看它是否能找到正常多项式的长度，比如下图所示的多项式：

图 11.4：一个复杂的多项式

这是一个 5 次多项式，意味着x的最高指数是5。曲线的方程如下：

图 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 来解决这个问题，方法是输入length of curve y = ... from –2 to 1，然后检查它是否是一个很好的近似值。但是使用 Python，有一种更直接的方法来计算曲线的长度，它不会遇到我们在导数中遇到的问题。事实上，它甚至不使用导数或积分。你可以简单地使用毕达哥拉斯定理找到曲线的微小部分的长度，并累加所有这些微小部分，如下图所示。我们知道宽度，我们对微小直角三角形的斜边感兴趣。我们可以计算高度，即x处的函数和x处的函数之间的差，再加上宽度：

图 11.6：找到曲线的微小部分的长度

在前面图表中显示的直角三角形的斜边如下：

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

我们所要做的就是遍历从a到b的区间，计算所有这些长度。以下是如何在 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值之间的曲线长度：

图 11.8：曲线方程

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

执行以下步骤：

1. 首先，我们需要使用前述方程创建一个circle函数：

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

注意

这又是一个半圆。这次，我们的curve_length2函数将毫无问题地累加弧线的微小切片。

1. 现在，我们将在该曲线上运行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

1. 现在，我们打印函数的输出，从x = -1到x = 1进行测量：

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

输出如下：

3.1415663562164773

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

注意

要访问此特定部分的源代码，请参阅packt.live/3gkI5Qi。

你也可以在packt.live/3eVpSbz上在线运行这个例子。

练习 11.02：找到正弦波的长度

数学和科学中一个非常重要和有用的函数是正弦波。它在 0 和 2π之间完成一个周期，如下图所示：

图 11.9：正弦波的一个周期

测量它的波长（2π）和振幅（它上下移动的距离，即 1 个单位）很容易，但实际曲线有多长呢？在这个练习中，我们将找到从 0 到 2π的正弦波的长度。

执行以下步骤：

1. 我们将再次使用我们的curve_length2函数，但现在我们必须从math模块导入我们的sin和pi函数：

from math import sin, pi

1. 我们已经编写了curve_length2函数，它将累加曲线的所有部分。我们只需要告诉它要使用的函数，以及开始和结束的x值：

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

输出如下：

7.640391636335927

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

注意

要访问此特定部分的源代码，请参阅packt.live/3dUy3nk。

您还可以在packt.live/2VFy2xd上在线运行此示例。

螺旋的长度

那么，极坐标中的螺旋怎么样，其中r，即与原点的距离，是与与x轴成角度（θ）的函数？我们不能使用我们的x和y函数来测量以下图中显示的螺旋：

图 11.10：阿基米德螺旋

在前面的图中，我们有一个从（5,0）开始，绕中心旋转 7.5 圈，最终到达（11,π）的螺旋。该曲线的公式为r(θ) = 5 + 0.12892θ。旋转的弧度数是 7.5 乘以 2π，即 15π。我们将使用与前一节相同的思路：我们将找到从r(θ)到r(θ+step)的直线长度，step 是中心角的一个小步长，如下图所示：

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

三角形中央角的对边就像我们积分问题中的切片或前一个曲线长度程序中三角形的斜边一样。这次，它不是直角三角形，所以我们不能使用斜边。但是对于这个问题，有一个称为余弦定律的公式。在三角形 ABC 中，角 C 的对边长度如下：

图 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. 我们不知道这个螺旋的公式，但我们知道半径在 12 个旋转中增加了 13 个单位（从 3 到 16）。这意味着对于角度θ的每个 2π增加，半径增加 13/12 个单位。因此，我们将 13/12 除以 2π。半径的增加可以表示如下：

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

1. 我们可以用 Python 这样表示：

def r(theta):
    return 3 + 0.1724*theta

1. 我们可以检查以确保r(0) = 3和r(24π) = 16：

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

输出如下：

3.0 15.998653763493127

1. 现在，我们只需将其放入我们的螺旋函数中：

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

输出如下：

716.3778471288748

在这个练习中，我们仅通过知道曲线的起始值和结束值以及它围绕中心旋转的次数，就轻松找到了这个螺旋曲线的长度，即716.3778471288748。

注意

要访问此特定部分的源代码，请参阅packt.live/2YT70EH。

您还可以在packt.live/2YV4wFT上在线运行此示例。

练习 11.04：找到卷中的绝缘长度

您被要求找到所示卷中绝缘的（近似）长度：

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

您测量卷起并发现中心是一个直径为 4 英寸的空圆（因此r(0) = 2）。卷的外径为 26 英寸。您从中心到外部计算层数，并估计螺旋需要 23 个半转，因此r(2π23.5)= 26/2 = 13*。

执行以下步骤计算长度：

1. 使用前面的数据计算方程：

图 11.15：计算半径的公式

1. 这就是螺旋图的样子：

图 11.16：绝缘卷的图表

1. 现在，将我们的r代码更改为这个螺旋并不难：

def r(theta):
    return 2 + 0.0745*theta

1. 现在，我们可以在这个函数上运行我们的螺旋函数，从0到2π23.5：

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

输出如下：

1107.502879450013

1,107 英寸只是超过 92 英尺的绝缘。

注意

要访问此特定部分的源代码，请参阅packt.live/2VE9YKZ。

您还可以在packt.live/31D43tG上在线运行此示例。

练习 11.05：找到阿基米德螺旋的长度

对于这个练习，您已经得到了阿基米德螺旋的方程。找到从θ=0到θ=2π的螺旋长度：

图 11.17：阿基米德螺旋的方程

注意

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

执行以下步骤找到长度：

1. 我们只需用指数函数重新定义r：

from math import e
def r(theta):
    return 2*e**(0.315*theta)

1. 然后，我们从0到2π运行螺旋函数：

spiral(r,0,2*pi)

输出如下：

41.518256747758976

这个螺旋的长度是41.518256747758976。

注意

要访问此特定部分的源代码，请参阅packt.live/2VEtjfo。

您还可以在packt.live/2VHasQN上在线运行此示例。

表面积

让我们学习如何将这个从二维转换为三维，并计算 3D 表面的面积。在第十章，使用 Python 进行基础微积分中，我们学习了如何计算旋转表面的面积，但这是一个第三维度z是x和y的值的函数的表面。

公式

解决这个问题的传统代数方法是通过对表面进行双重积分：

图 11.18：计算表面积的公式

这里，z = f(x,y)或(x,y,f(x,y))。那些花括号 d 是增量，意味着我们将处理偏导数。偏导数是导数，但只针对一个变量，即使函数依赖于多个变量。这是一个返回函数f在特定点(v,w)处相对于变量u的偏导数的函数。根据我们感兴趣的变量是x还是y，函数将朝着那个方向迈出微小的一步并计算导数，就像我们已经做过的那样：

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

在代码中有一个try...except块，以防抛出ValueError。如果斜率太大，就像垂直线一样，就会发生这种情况。如果发生这种情况，它将忽略它并继续进行。

现在，我们需要一个 3D 向量和一个cross函数来计算面积公式中的叉积。叉积给出垂直于给定向量的向量的长度，但也给出由给定向量形成的平行四边形的面积。

图 11.19：两个向量的叉积

如果您知道向量之间的角度，可以使用它来找到叉积：

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

如果不这样做，就像我们的情况一样，可以使用 3D 向量来表示向量在每个方向上的位移，x，y和z。例如，假设我们有两个向量，u = 2i + 3j + 4k和v = 5i + 6j + 7k。它们由它们在三个维度中的位移定义。i部分是x方向的位移，j部分是y方向的位移，k部分是z方向的位移。好消息是，会有一些零来简化事情。要交叉两个向量，我们可以将它们放入矩阵中，并对它们进行操作，如下矩阵的行列式：

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

我们将编写一个函数来执行两个 3D 向量的操作。我们需要放入的只是i，j和k的系数。因此，如果u = ai + bj + ck和v = di + ej + fk，我们将得到以下结果：

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

让我们使用列表来表示向量，所以u = [a,b,c]和u[0] = a，u[1] = b和u[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]]

我们编写了一个长的文档字符串，以清楚地说明函数的用途，如何放入值以及我们将得到什么输出。让我们检查一下，以确保我们得到正确的输出：

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

输出如下：

[-3, 6, -3]

可以了。现在，我们需要编写一个函数来找到 3D 向量的大小，因为这将给我们平行四边形的面积。这只是将勾股定理扩展到三维。因此，如果u = ai + bj + ck，则向量u的大小是：

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

这是半圆的样子，其表面由平行四边形近似。更多的平行四边形应该意味着更准确的近似：

图 11.23：使用更多的平行四边形

我们的面积函数将循环遍历网格中的所有x和y点，计算每个点的偏导数，并使用叉积来找到该点处平行四边形的面积：

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。然后，我们计算dx和dy，x和y的微小变化，因为我们将表面分成相等的切片。try...except块只是忽略(pass)当切线的斜率是垂直的时候，会出现无限的偏导数的错误，就像我们在图 11.3中看到的那样。如果没有错误，它将添加在该点由偏导数形成的平行四边形的面积。现在，我们使用 1,000 个点在每个方向上运行半球的面积函数，并得到一个相当准确的近似值。我们知道半径为 1 的球体的表面积的一半是 2π，或 6.28：

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

输出如下：

Area of hemisphere: 6.210356913122

现在，让我们快速进行一个基于这个概念的练习。

练习 11.06：寻找 3D 表面的面积-第 1 部分

现在，让我们找到一个复杂表面的面积，这将很难使用代数方法找到。考虑以下方程的表面：

图 11.24：表面的方程

表面显示在以下图像中：

图 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))

1. 运行程序以查看输出：

Area of wave surface: 608\. 2832236305994

从前面的代码中，我们可以清楚地看到，使用 Python 几行代码就可以轻松找到甚至复杂表面的面积有多容易。

注意

要访问此特定部分的源代码，请参阅packt.live/3gwd6kr。

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

练习 11.07：找到 3D 曲面的面积-第 2 部分

找到曲面的面积。

这就是曲面的样子：

图 11.26：另一个 3D 曲面

执行以下步骤找到面积：

1. 定义我们的曲面函数以返回表达式：

def surface(x,y):
    return 3*cos(x)+2*cos(x)*cos(y) 

1. 运行surface函数以获得值：

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

输出如下：

Area of surface: 99.80676808568984

这个 3D 曲面的面积是99.80676808568984。

注意

要访问此特定部分的源代码，请参阅packt.live/2VCaObq。

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

练习 11.08：找到曲面的面积-第 3 部分

找到曲面的面积。

这就是曲面的样子：

图 11.27：曲面

执行以下步骤找到面积：

1. 定义我们的曲面函数以返回新表达式：

def surface(x,y):
    return sqrt(1+sin(x)*cos(y))

1. 运行surface函数：

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

输出如下：

Area of surface: 42.80527549685105

这个曲面的面积是42.80527549685105。

注意

要访问此特定部分的源代码，请参阅packt.live/3gwdLlV。

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

无限级数

数学家经常遇到对他们来说太复杂以至于无法解决或处理的函数，近似一直是数学中的重要组成部分。对于试图代数地进行导数和积分的数学家来说，许多表达式没有漂亮的解、导数、积分等。一般来说，科学家在现实生活中遇到的微分方程没有代数解，因此他们必须使用其他方法。稍后会详细介绍微分方程，但有一类重要的近似使用简单函数来近似难函数。

多项式函数

解决多项式方程很容易，可以求导和积分，比如y = x2，甚至以下方程：

图 11.28：多项式方程

所有项都是一个接一个地添加（或相减），没有三角、对数或指数函数参与，使事情变得困难。以下是用简单多项式近似难函数的公式：

图 11.29：泰勒级数

这个公式被称为泰勒级数：任何函数（可导）都可以使用一定精度在特定点用一定页数的多项式来近似。

级数

数学家有一种表示将一堆遵循某种模式的数字相加的符号：

图 11.30：计算级数的公式

看起来像一个E的大符号实际上是希腊字母 sigma，或者S，表示数字的总和。sigma 下面的方程是变量从哪里开始（在这种情况下是 1），sigma 上面是i的最后一个整数值（在这种情况下是 10）。在 sigma 右边是对变量的操作的表达式。在这种情况下，我们只是将i，变量，从 1 加到 10。这几乎就是您在 Python 中编写列表推导式的方式。如下：

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

列表推导式中的第一个项是在 sigma 级数表达式中看到的内容-在这种情况下是i。例如，直到n的整数平方和的级数如下：

图 11.31：1 到 n 的整数平方和的级数

在 Python 中，我们会这样写：

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

一个古老但有用的级数是 arctangent 级数。它计算具有给定正切的角度（以弧度表示）；例如：

图 11.32：arctangent 级数的方程

根据前述方程，arctan 的方程如下：

图 11.33：arctan 的方程

该系列按照以下模式计算：

图 11.34：arctan 级数的方程

这是 sigma 表达式：

图 11.35：一个 sigma 表达式

通过插入x的正切，我们可以计算出角度的近似值：

图 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

这非常接近。

收敛

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

图 11.37：正切的三角函数

根据前述方程，arctan 的方程如下：

图 11.38：计算 arctan 的公式

他们认为将x替换为 arctan 级数中的1会使计算π变得轻而易举。以下是前几个项：

图 11.39：将 x = 1 代入 arctan 级数

这个表达式给出了π的近似值：

图 11.40：寻找π的近似值的方程

我们只需编写 sigma 右侧的代码，添加n的范围代码，并将它们相加：

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

我们可以在输出中显示逼近π的进展：

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 年，英国数学家和天文学家约翰·马钦使用他改进的级数计算了π的 100 位小数。以下是该级数：

图 11.41：一个 arctan 函数

使用前述 arctan 函数来计算π的 10 个正确数字。按照以下步骤进行：

1. 只需调用我们的 arctan 函数。10 个项应该足够了：

print(4*(4*arctan(1/5,10)-arctan(1/239,10)))

1. 运行上述代码以查看输出：

3.1415926535897922

使用 10 个项可以得到一个相当不错的近似值。它甚至给出了超过 10 个正确的数字。

注意

要访问此特定部分的源代码，请参阅packt.live/3dPjVvD。

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

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

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

图 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)])))

1. 它收敛得更快吗？运行前面的代码以查看输出：

3.1320765318091053
3.1406380562059946
3.1415916986605086

看起来它似乎并没有更快地收敛。100 万项后，您仍然只有五位正确的小数。

注意

要访问此特定部分的源代码，请参阅packt.live/2NRnnLD。

您还可以在packt.live/38lHXgm上在线运行此示例。

20 世纪的公式

这是卓越的印度数学家拉马努金用来近似π的公式：

图 11.43：拉马努金近似π的公式

以下是如何在 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 部分

确定以下幂级数的收敛区间：

图 11.44：幂级数

执行以下步骤：

1. 将总和输入 Python：

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

由于我们不能使用数字“无穷大”，我们找到了从 n = 1 到 100 万的所有项的和。

1. 运行所有从-10 到 10 的整数，看看是否收敛到一个数字：

for x in range(-10,11):
    print(x,mystery_sum(x))

1. 当您运行此代码时，将会出现OverflowError：

OverflowError: int too large to convert to float

这意味着一些数字变得很大，这正是我们预期的。我们需要添加一个条件，以便如果出现错误，它将简单地返回Infinity。这是通过try...except块完成的。

1. 让我们告诉 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"

1. 现在，输出给我们一些无穷大和一些实际值：

-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。这意味着如果x在该区间内，我们可以使用该级数获得有用的值。否则，我们无法使用它。

注意

要访问此特定部分的源代码，请参阅packt.live/38k30A2。

您还可以在packt.live/31AtmMU上在线运行此示例。

练习 11.12：确定收敛区间-第 2 部分

确定以下幂级数的收敛区间：

图 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))

1. 以下是一些输出：

-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

在x在-25 和 25 之间的所有输出都保持很小（在 0 和 600 之间），无论我们使用多少项，因此我们将称之为收敛区间-25 < x < 25。

注意

要访问此特定部分的源代码，请参阅packt.live/38pwuwC。

您还可以在packt.live/2YS46jl上在线运行此示例。

练习 11.13：找到常数

在这个练习中，我们将在 Python 中表示一个无限级数并找到总和。我们将使用一个著名的常数，它被定义为该级数的和：

图 11.46：级数的总和

这个著名常数的值是多少？让我们按照以下步骤来确定这个值：

1. 导入阶乘模块并将前述方程转换为 Python，如下所示：

from math import factorial
print(sum([1/factorial(n) for n in range(10000)]))

1. 运行前面的代码以查看输出：

2.7182818284590455

著名的常数是e，自然对数的底数。

注意

要访问此特定部分的源代码，请参阅packt.live/2AoyubH。

您还可以在packt.live/2BZ4aVw上在线运行此示例。

活动 11.01：寻找表面的最小值

机器学习中的一个重要任务是最小化函数。当您训练神经网络时，您正在改变矩阵或张量中的值，以查看哪些值提供了更好地逼近您的测试数据。在网络的每个值处，您都可以看到它对您的错误有多大贡献。这听起来像是表面上不同点的偏导数，不是吗？

这个例子是梯度下降的过程。让我们考虑一下，我们想要找到我们函数的最小值的位置。我们表面上的每个点都有一个偏导数，我们可以使用这些偏导数向更低的值移动一点。我们将从一个随机的地方开始，计算该点的偏导数，然后沿着降低z值的方向移动，也就是上下值。因此，如果z关于x的偏导数（我们称为dz_dx）是负的，这意味着z随着x的增加而减少，我们将希望向正x方向移动。如果dz_dx是正的，这意味着z随着x的增加而增加，所以我们将希望朝相反的方向移动，因此我们将向负x方向移动。我们将对y方向做同样的事情。这将如下所示：

图 11.47：下降到最小值的路径

这个活动的第一部分是创建一个找到表面最小点的函数。可以通过以下步骤编写此函数：

1. 编写一个函数，该函数将在表面上创建一个随机的(x, y)位置。您可以调用random模块的uniform函数来生成这些值。

2. 计算z关于x和y的偏导数。

3. 通过偏导数的负值乘以一个微小的step量来改变x和y，以防偏导数很大。

4. 计算这个新位置的偏导数，并保持循环，直到偏导数都非常小（小于 0.0001），或者位置偏离表面。

5. 在一堆随机位置上运行函数，将最小z值保存到mins列表中。

6. 最后，打印mins列表的最小值。

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

1. 在表面上测试您的函数。您的函数应该发现最小值为 0，在点(0, 0)处。

2. 一旦您对您的函数有信心，就可以使用它来确定的最小值，其中-1 < x < 5和-1 < y < 5。

您会发现，根据您的起始点，您的函数将收敛到不同的最小点 - 局部最小值和全局最小值。

注意

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

总结

在上一章中，我们学习了导数和积分的力量，因此在本章中，我们基于这些工具来解决一些非常困难的问题，比如螺旋线的长度和三维表面的面积。我们甚至通过引入偏导数将导数和积分扩展到三维空间。在微积分课上，我们需要使用大量的代数来使用这些工具，但是通过使用 Python，我们对情况进行了建模并测试了我们的函数。我们创建了包含我们变化数值的变量，并在必要时循环计算数百万次。对于以前世纪的数学家来说，这似乎就像是某种魔法灯。

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

FAB62

RUC47

第十二章：12.使用 Python 进行中级微积分

概述

在本章结束时，你将能够解决涉及变量变化的方程的问题。在本章中，你将使用数值方法来模拟人口和温度，并使用微分方程来计算它们的过去值或预测它们的未来值。当你知道数字在特定范围内时，你将学会使用二分搜索来猜测和检查，以获得非常准确的解决方案。你还将模拟物体移动的情况，并在给定它们的速度微分方程时解决它们的未来位置。

介绍

数学学生经常抱怨代数和几何中出现的问题没有现实世界的应用，比如因式分解多项式或者角的二等分，但对于微分方程却不能这样说。使用本章中将学到的工具，你将能够用微分方程模拟和解决物理、电子和工程中的真实问题。Python 是数学家和科学家的完美工具，他们想要进行数字计算和解决问题，但又不想为此再获得计算机科学学位。Python 是最受欢迎的编程语言之一，因为它易于使用且没有不必要的抽象。

到了 17 世纪，数学家们用数学方程模拟了物体的运动，并把目光投向了外层空间的行星。牛顿模拟了它们的运动，他提出的方程不仅涉及未知数，还涉及这些数字的变化。例如，他的方程不仅包含一个未知角度，还包含这个角度的变化（它的角速度），甚至是角度变化的变化（它的加速度）！当时没有工具来解决这些方程，所以他不得不自己发明这些工具。这些工具后来被称为微积分。

微分方程

在数学课上解方程通常涉及一个未知数x。方程隐藏了这个值，但给了你一些提示，告诉你如何找到这个值，比如。但是要解微分方程，你只能得到关于函数的导数的信息，然后期望你找到这个函数。可能是像下面这样简单的东西：

图 12.1：找到导数为 2 的函数

这意味着找到一个导数为 2 的函数。这也可以写成如下形式：

图 12.2：表示函数导数的另一种方式

通过简单的积分，我们可以找到适用于这个方程的函数，因为我们知道函数y = 2x的导数是 2。事实上，许多相关函数，比如y = 2x + 1，y = 2x + 2，y = 2x + 3等等，它们的导数都是 2。所以，我们写出一个一般形式，即y = 2x + C。

当我们没有太多线索时，事情就变得更加复杂，就像在这个方程中：

图 12.3：函数值为函数本身的导数

这是要求一个导数为函数本身的函数。

为了理解微分方程是如何使用的，让我们从简单的函数开始，以及涉及到现实世界中的东西，比如金钱。

利息计算

微分方程研究中有一个关键工具起源于中世纪的利息计算研究。让我们来看下面的练习。

练习 12.01：计算利息

一个储蓄账户每年支付 2%的利息。如果投资了 3500 美元，5 年后账户中有多少钱？

简单利息的公式如下：

图 12.4：简单利息公式

在这里，I是利息，P是本金或原始投资金额，r是利率或增长率，t是投资金额已投资的时间。根据这个公式，金额的利息为I = (3500)(0.02)(5) = 350 美元。按照以下步骤完成此练习：

1. 这是一个很好的机会来开始一个程序，它将接受初始金额、利率和时间，并使用前面的公式输出利息收入：

def amount(p0,rate,t):
    """Returns the amount after t
    years starting at p0 and growing
    at the given rate per year"""
    return p0*rate*t

1. 如您在amount函数的文档字符串中所看到的，它将接受一个起始金额和增长率，并返回给定年数后的投资金额。让我们看看 1-5 年内的利息收入：

for i in range(1,6):
    print(i,amount(3500,0.02,i))

以下是输出：

1 70.0
2 140.0
3 210.0
4 280.0
5 350.0

但这并不是利息的真正工作方式。每年几次，我们计算该年份的利息收入，将其加到本金中，新的本金更高。下一次的利息计算是在更高的数字上，因此称为复利。给定每年 n 次复利的t年后的金额的公式如下：

图 12.5：计算 t 年后的金额公式

1. 让我们将amount函数更改为以下内容：

def amount(p0,rate,t,comps):
    """Returns the amount after t
    years starting at p0 and growing
    at the given rate per year
    compounded comps times per year"""
    for i in range(int(t*comps)):
        p0 += p0*rate/comps
    return p0

在这个函数中，我们添加了按复利次数给出的年份的利息收入。如果我们每年只计算一次复利，看起来是这样的：

for i in range(1,6):
    print(i,amount(3500,0.02,i,1))

这就是我们得到的：

1 3570.0
2 3641.4
3 3714.228
4 3788.51256
5 3864.2828112

因此，在 5 年结束时，我们赚了 364 美元，而不仅仅是 350 美元。即使利率相同，复利更频繁也会使金额增长更快。如果我们将复利更改为每年 12 次（每月复利），我们将在 5 年后得到 3867 美元，比年复利多一点。

注意

要访问此特定部分的源代码，请参阅packt.live/3dUWz7C。

您也可以在packt.live/3iqUKCO上在线运行此示例。

练习 12.02：计算复利-第 1 部分

在一个年利率为 5.5%的储蓄账户中投资了 2000 美元，按月复利。要将金额增长到 8000 美元需要多长时间？按照以下步骤来解决这个问题：

1. 我们将使用我们从上一个练习中的amount函数打印出投资的前 5 年：

for i in range(1,6):
    print(i,amount(2000,0.055,i,12))

输出如下：

1 2112.815720771071
2 2231.9951349686903
3 2357.8972049231984
4 2490.9011412619493
5 2631.4075450724245

1. 5 年后，金额只有 2631 美元。要达到 8000 美元，我们必须走 20 或 30 年：

for i in [5,10,15,20,25,30]:
    print(i,amount(2000,0.055,i,12))

输出如下：

5 2631.4075450724245
10 3462.1528341320413
15 4555.167544964467
20 5993.251123444263
25 7885.343112872511
30 10374.775681348801

在 25 到 30 年之间的某个时候，我们将达到 8000 美元。更精确的方法是更聪明地猜测。

1. 我们将范围减半，并根据我们得到的结果猜测更高或更低。例如，25 年和 30 年的平均值是 27.5，因此我们输入以下内容：

print(amount(2000,0.055,27.5,12))

以下是输出：

9044.814313545687

因此，我们将在 27.5 年内达到 9000 美元。达到 8000 美元的时间必须少于这个时间。

1. 我们将计算 25 和 27.5 的平均值并将其代入：

def average(a,b):
    return (a+b)/2
print(amount(2000,0.055,average(25,27.5),12))

以下是输出：

8445.203624219383

1. 让我们编写一个程序，直到找到答案为止。这称为二分搜索。让我们创建一个bin_search函数，它将使用我们正在使用的函数的名称，我们正在搜索的范围的下限和上限以及目标输出（在本例中为 8000）作为参数：

def bin_search(f,lower,upper,target):
    for i in range(20):
        avg = average(lower,upper)

1. 这是关键的一行。它将平均值插入函数中，使用所有其他必需的参数，并将输出分配给guess变量。我们将检查该变量是否等于我们的目标，或者我们是否需要猜测更高或更低：

        guess = f(2000,0.055,avg,12)
        if guess == target:
            return guess
        if guess > target:
            upper = avg
        else:
            lower = avg
    return avg

1. 我们将我们范围的下限和上限以及我们的目标数字代入我们的函数，以获得我们的近似值：

bin_search(amount,25,30,8000)

输出如下：

25.333333015441895

1. 看起来我们将在25 年零 4 个月内达到 8000 美元。让我们检查一下：

amount(2000,0.055,25.334,12)

果然，复利后的余额超过了 8000 美元：

8030.904658737448

我们将再次使用二分搜索，但现在，让我们使用我们的代码来找到一个在微分方程中经常出现的相当重要的数学常数。

注意

要访问此特定部分的源代码，请参阅packt.live/3iq95PV

您还可以在packt.live/2BpdbHI上在线运行此示例

练习 12.03：计算复利-第 2 部分

如果您以 100%的利率投资$1，连续复利 1 年，您会赚多少？

请记住，复利的频率越高，最终金额就越高。您认为会是多少？$1.50？$2？本金、利率和时间都是 1，但comps变量是什么？按照以下步骤完成这个练习：

1. 为了近似连续复利，我们将每秒复利一次（每年365246060*次）：

print(amount(1,1,1,365*24*60*60))

输出如下：

2.7182817853606362

大约是$2.72。这个数字，2.71828…，是自然对数的底数e。它对于模拟自然界中的人口非常有用，因为动物、植物和微生物不会等到月底才繁殖-它们是持续不断地繁殖。因此，当利息连续复利或人口自然增长时，我们将使用这个公式：

图 12.6：计算复利的公式

1. 让我们创建一个函数来快速完成这个任务。首先，我们需要从math模块中导入e以进行连续复利：

from math import e

1. 创建一个pert函数，它将插入初始金额或人口、增长率和时间，并返回最终金额：

def pert(P,r,t):
    return P*e**(r*t)

我们将在本章中多次使用这个函数。现在，让我们回答更多的投资问题。

注意

要访问此特定部分的源代码，请参阅packt.live/2D2Q1r0

您还可以在packt.live/31G5pDQ上在线运行此示例

练习 12.04：计算复利-第 3 部分

一个人以每月复利的 18%年利率借了$5,000。1 年后这个人会欠多少钱？按照以下步骤完成这个练习：

1. 我们可以将其放入我们的函数调用中：

amount(5000,0.18,1,12)

输出如下：

5978.090857307678

为了比较，让我们看看如果利息是连续复利会发生什么。

1. 我们将使用我们的pert函数输入P = 5000，r = 0.18和t = 1作为值：

print("Continuous:",pert(5000,0.18,1))

得到的金额如下：

5986.096815609051

注意

要访问此特定部分的源代码，请参阅packt.live/31ES5Qi

您还可以在packt.live/3f5j0s4上在线运行此示例

练习 12.05：成为百万富翁

如果您以每日复利 8%的年利率投资$1,000，要成为百万富翁需要多长时间？如果初始金额是$10,000 呢？按照以下步骤完成这个练习：

1. 首先，让我们定义bin_search函数，如下所示：

def bin_search(f,lower,upper,target):
    for i in range(20):
        avg = average(lower,upper)
        #Be sure to change this line
        #if the principal, rate or
        #compounding changes:
        guess = f(1000,0.08,avg,365)
        if guess == target:
            return guess
        if guess > target:
            upper = avg
        else:
            lower = avg
    return avg

1. 让我们猜一些野生的猜测，看看如果$1,000 投资这些年，我们会得到多少：

for i in [10,20,30,40,50]:
    print(i,amount(1000,0.08,i,365))

这是输出：

10 2225.34584963113
20 4952.164150470476
30 11020.277938941583
40 24523.929773205105
50 54574.22533744746

1. 50 年后，您仍然只有$54,000，而不是一百万。但是 100 年后，您将拥有近 300 万：

amount(1000,0.08,100,365)

这是输出：

2978346.0711824815

1. 答案必须在 50 和 100 之间。看起来是我们二分搜索的任务：

print(bin_search(amount,50,100,1000000))

我们得到以下输出：

86.3588809967041

1. 这表明在 86.36 年后，我们将拥有 100 万美元。如果初始投资是$10,000，那么在bin_search函数中更新guess变量：

        guess = f(10000,0.08,avg,365)

这是我们将打印所需输出的方法：

for i in [10,15,20,30,40,50,60]:
    print(i,amount(10000,0.08,i,365))

输出如下：

10 22253.458496311334
15 33196.803971077774
20 49521.64150470513
30 110202.77938941623
40 245239.2977320514
50 545742.2533744735
60 1214465.2585152255

1. 因此，我们在 50 到 60 年之间就能达到 100 万美元。让我们在我们的二分搜索函数中将1000改为10000并检查一下：

print(bin_search(amount,50,60,1000000))

我们得到以下输出：

57.57260322570801

超过 57.57 年才能达到一百万美元。

因此，我们开始学习微分方程是通过研究复利来学习的。一定金额的钱每年、每月或每天都会有利息率。初始金额

注意

要访问此特定部分的源代码，请参阅packt.live/31ycoPg

您也可以在packt.live/2NMT9sX上在线运行此示例

现在，我们将同样的推理扩展到人、动物、细菌和热量的数量，这些数量不断变化，或者说是连续的。

人口增长

微分方程对于找到特定时间的人口、动物和细菌的数量的公式非常有用；例如：

图 12.7：计算时间 t 的微分方程

这个微分方程意味着y的变化速率与y成比例，或者说人口的增长与其数量成比例。这就是人口增长率的定义：人口的一部分或百分比。解决方案类似于涉及连续复利的利息问题：

图 12.8：计算变化率的微分方程

练习 12.06：计算人口增长率-第 1 部分

在 1980 年代，肯尼亚的年人口增长率为 4%。以这个速度，人口翻倍需要多长时间？按照以下步骤完成这个练习：

1. 无论初始人口是多少，我们都在寻找使因子ert 等于 2 的t。我们可以使用我们的pert函数和二项式搜索函数，稍作调整：

def bin_search(f,lower,upper,target):
    for i in range(40):
        avg = average(lower,upper)
        guess = f(1,0.04,avg)
        if guess == target:
            return guess
        if guess > target:
            upper = avg
        else:
            lower = avg
    return avg

1. 我们正在寻找时间t，以 4%的增长率将我们的初始人口从 1 增加到 2。我们估计这个时间会在 1 到 100 年之间：

print(bin_search(pert,1,100,2))

输出如下：

17.32867951408025

我们可以用代数来验证这一点。我们取方程两边的对数并解出t：

图 12.9：解决时间(t)的方程

1. 这意味着在 17 年多一点的时间内，肯尼亚的人口将翻倍。我们可以用我们的amount函数来验证这一点。1989 年，肯尼亚的人口为 2100 万：

print(amount(21000000,0.04,17.3,1000000))

以下是输出：

41951845.46179989

是的，每年使用一百万次复利，人口在 17.3 年内增长到了将近 4200 万。

作为对此的回应，肯尼亚政府大力推广计划生育。这有效果吗？

注意

要访问此特定部分的源代码，请参阅packt.live/2BxsfCT

您也可以在packt.live/2Zuoy9c上在线运行此示例

练习 12.07：计算人口增长率-第 2 部分

2010 年，肯尼亚的人口为 4200 万。到 2019 年，人口为 5250 万。该范围的年人口增长率是多少？

我们可以再次使用二分搜索函数，返回给定初始人口（以百万计）、时间t和 9 年后的目标人口（以百万计）的增长因子r。

在bin_search函数中，将时间更改为9：

        guess = f(1,avg,9)

然后，我们将找到这 9 年的年增长率。我们知道它在 0 到 2 之间：

print(bin_search(pert,0,2,52.5/42))

打印出的值如下：

0.024793727925498388

计划生育项目一定是有效的！肯尼亚将其人口增长率降低到了 2.5%。

注意

要访问此特定部分的源代码，请参阅packt.live/3eWKzDW

您也可以在packt.live/31EKPUq上在线运行此示例

放射性材料的半衰期

与人口问题类似，半衰期问题涉及一个群体，但其中一半是放射性材料的原子，其中一半的原子随着时间变成了不同物质的原子。例如，碳-14 衰变成氮-14，大约需要 5730 年才能使一半的碳衰变。这使得放射性碳测年成为从考古学到检测伪造艺术品的关键工具。

练习 12.08：测量放射性衰变

镭-226 的半衰期为 1600 年。在给定样本中，800 年内将消失多少镭？

意思是“物质的衰变速率与物质的数量成比例”的微分方程表达如下：

图 12.10：用于计算物质衰变速率的微分方程

解决方案与我们的人口问题类似，只是衰减因子是负的，因为数量减少：

图 12.11：使用负衰变因子计算变化率

这意味着最终金额等于时间的初始金额，e，乘以衰减因子，r，和时间，t的乘积。我们可以像解决人口问题一样使用我们的二分搜索函数。我们正在寻找在 1,600 年内使我们的人口减半的增长率r。按照以下步骤完成这个练习：

1. 在bin_search函数中的guess =行中将t更改为1600：

        guess = f(1,avg,1600)

1. 然后，搜索增长因子，我们认为它将在-2 和 0 之间。我们的目标金额是起始金额的一半：

print(bin_search(pert,-2,0,0.5))

以下是输出：

-0.0004332169864937896

1. 这就是镭-226 的衰变因子r。我们要做的就是将其插入我们的pert函数中，以找出 800 年后剩下的样本的百分比：

pert(1,-0.0004332,800)

以下是输出：

0.7071163910309745

因此，大约 71%的样本在 800 年后仍然存在。

注意

要访问此特定部分的源代码，请参阅packt.live/2YSzQ84

您还可以在packt.live/2ByUwJj上在线运行此示例

练习 12.09：测量历史文物的年龄

对一块布进行了放射性碳测年。这意味着科学家们测量了有多少碳-14（半衰期 5730 年）衰变成了更稳定的同位素。他们发现剩下的碳-14 的数量是碳-13 的 10 倍。这块布有多大年龄？

如果碳-14 需要 5730 年使其数量减半，我们需要找到我们的 Pert 公式的速率r：

图 12.12：Pert 公式

按照以下步骤完成这个练习：

1. 我们使用我们的二分搜索函数来解决r：

def bin_search(f,lower,upper,target):
    for i in range(40):
        avg = average(lower,upper)

1. 这是更改的那一行。如果我们在pert函数中放入一个起始金额为1，r将是avg，5730将是目标时间：

        guess = f(1,avg,5730)
        if guess == target:
            return guess
        if guess > target:
            upper = avg
        else:
            lower = avg
    return avg
print(bin_search(pert,-2,0,0.5))

以下是输出：

-0.00012096809405193198

r = -0.000120968，因此我们的 Pert 公式变为如下：

图 12.13：在 Pert 公式中替换 r 的新值

这意味着x克的碳-14 衰变了，剩下的是 10x 克，是整个样本的 10 倍。因此，衰变量是整个样本的 1/11 或 0.091。结束金额是 1-0.091。这使我们的 Pert 方程如下：

图 12.14：带有结束金额的 Pert 方程

1. 我们方程中唯一未知的是t，因此我们正在更改我们的bin_search函数，以便有策略地猜测和检查正确的t。返回到您的bin_search函数；开头应该是这样的：

def bin_search(f,lower,upper,target):
    for i in range(40):
        avg = average(lower,upper)

1. 这是我们正在更改的行。我们将 1 代入原始量，长小数是我们的r，时间范围的平均值用于时间。目标是样本的 0.091，这将保持猜测和平均值，直到返回确切的年数以达到目标值：

        guess = f(1,-0.000120968,avg)
        if guess == target:
            return guess

1. 由于它是一个递减函数，如果猜测小于目标值，我们将会超过目标值，upper数字将被替换为平均值：

        if guess < target:
            upper = avg
        else:
            lower = avg
    return avg
print(bin_search(pert,1,100000,0.91))
print(pert(1,-0.000120968,5730))

1. 请注意，我们更改了if guess < target:行。我们正在寻找从 1 到 0.91 的衰减量所需的年数，以给定的速率。我们怀疑它在 1 到 100,000 年之间。第二个print行只是检查我们的pert函数确认在 5,730 年后，剩余量正好是原始量的一半。当我们运行代码时，这是输出：

779.633287019019
0.5000002702800457

根据我们的计算，这块布大约有780 年的历史。

因此，我们最初编写此代码是为了测量投资中剩余的金额，该投资以给定的利率增长一段时间。在本节中，我们将此应用于物体中放射性物质的剩余量，该物质以已知速率衰变，时间未知。这就是科学家计算考古文物年龄的方法。

注意

要访问此特定部分的源代码，请参阅packt.live/3eOJJJv。

您也可以在packt.live/38mESgn上在线运行此示例。

接下来，我们将使用相同的思路，但将其应用于物体的温度变化，比如一杯咖啡或人体的温度。

牛顿冷却定律

你是否曾经想过在警察节目中带着乳胶手套的犯罪现场调查员（CSI）如何判断受害者的死亡时间？艾萨克·牛顿被认为是发现物质冷却遵循微分方程的人：

图 12.15：温度变化速率的微分方程

看看这个微分方程与我们以前看到的微分方程略有不同？这不是物质温度变化速率与物质温度成比例，而是说“物质温度变化速率与物质温度与环境温度之间的差值成比例”。因此，如果一杯热咖啡放在热的房间里，它的温度变化速度会比放在非常冷的房间里慢。同样，我们知道警察节目中受害者的体温起始温度：98.6°F。

练习 12.10：计算死亡时间

一名调查员到达犯罪现场并测量环境和尸体的温度。如果环境温度为 65°，尸体温度为 80°，调查员记录时间并等待一小时。尸体温度与环境温度的差为 15 度。一个小时后，环境温度仍为 65°，尸体进一步冷却至 75°。温度差现在为 10 度。受害者是何时死亡的？

有了这些信息，她可以建立以下方程：

图 12.16：计算死亡时间的方程

按照以下步骤完成这个练习：

1. 我们可以使用二分搜索来找出温度的衰减率。我们需要导入e并确保我们有pert和average函数：

from math import e
def pert(P,r,t):
    return P*e**(r*t)
def average(a,b):
    return (a+b)/2

1. 我们的bin_search函数的第一部分与以前相同：

def bin_search(f,lower,upper,target):
    for i in range(40):
        avg = average(lower,upper)

1. 这里的重要变化是：我们的原始量（温度差）为 15 度，我们想知道r，即我们 Pert 公式中的变化率：

        guess = f(15,avg,1)
        if guess == target:
            return guess
        if guess > target:
            upper = avg
        else:
            lower = avg
    return avg
print(bin_search(pert,-2,0,10))

这是输出：

-0.4054651081078191

这就是这种情况的衰减速率，所以我们知道了尸体温度与环境温度之间的初始差异（98.6-65），以及最终差异（10）和衰减速率。这是情况的图表：

图 12.17：冷却体的图表

我们只需要知道差异衰减到 10 所需的小时数。我们的方程如下：

图 12.18：差异衰减到 10 所需的小时数

1. 我们改变我们的二分搜索函数来获取时间：

def bin_search(f,lower,upper,target):
    for i in range(40):
        avg = average(lower,upper)
        guess = f(33.6,-.4055,avg)
        if guess == target:
            return guess
        if guess > target:
            upper = avg
        else:
            lower = avg
    return avg

但是，如果时间太短，差异将太大。绕过这个最简单的方法是将更高的t作为函数调用的较低部分，将较低的t作为搜索范围的上限。

1. 调查员推断时间必须在 0 到 5 小时之间：

print(bin_search(pert,5,0,10))

输出将如下：

2.9887570258370033

几乎正好 3 小时。这看起来接近前图中曲线的y值为 10 的时间。

1. 让我们在我们的pert函数中检查一下。从r = -0.4055和t = 3.0开始，差异为 33.6 度。希望最终得到 10：

pert(33.6,-0.4055,3)

以下是输出：

9.954513505592326

所以，现在，当明星侦探在凌晨 2:30 到达现场时，调查员可以说，“死亡时间大约是晚上 11:30。”

注意

要访问此特定部分的源代码，请参阅packt.live/38jN68K.

您还可以在packt.live/3gefegi.上在线运行此示例

练习 12.11：计算温度变化的速率

一杯温度为 175°F 的咖啡放在一个 72°F 的房间里。我们等待 15 分钟，测量咖啡的温度，发现它已经变成了 140°F。按照这个速度，从开始算起 1 小时后它的温度会是多少？按照以下步骤完成这个练习：

1. 差异从 103°（175-72）开始。在 0.25 小时内，它变为 68°（140-72）。现在，我们可以建立一个方程：

图 12.19：计算咖啡温度差异的方程

1. 我们可以改变我们的二分搜索函数以反映这种情况。将guess=行更改为bin_search函数中的以下内容：

        guess = f(103,avg,0.25)

1. 运行它，找出在-2 和 0 之间的r将给我们带来 68°的差异：

print(bin_search(pert,-2,0,68))

这是输出：

-1.6608851322143892

1. 太快了！将其放入我们的 Pert 公式中，P = 103和t=1：

pert(103,-1.6608851322143892,1)

以下是输出：

19.566987911888482

这是 1 小时的差异。如果房间温度为 72°，那意味着咖啡将是72 + 19.5 = 91.5°。

注意

访问此特定部分的源代码，请参阅packt.live/3gl5p0i.

您还可以在packt.live/2YTdCmw.上在线运行此示例

混合问题

在代数中，有一些需要计算的文字问题，您必须计算出需要添加多少材料到混合物中才能获得特定的浓度或数量。在微积分中，自然，问题必须更难：例如，混合物正在改变；材料正在进入混合物，材料正在流出。您必须找出在特定时间后有多少混合物或溶剂。让我们看看以下练习，以更好地理解这个概念。

练习 12.12：解决混合问题-第 1 部分

一个罐子里含有 82 加仑的卤水，其中溶解了 18 磅的盐。每分钟以 5 加仑的速度流入罐子的卤水中含有每加仑 3 磅的溶解盐。这种混合物通过搅拌保持均匀，以每分钟 2 加仑的速度流出罐子。39 分钟后罐子里有多少盐？

正如你可以想象的，这种问题会导致一些复杂的微分方程，只有经过几页的代数运算，你才能得到一个方程（通常涉及到* e*的某个幂），然后你可以将时间代入并得到最终的数量。然而，使用编程，我们只需从给定的起始溶液开始，并添加和减去问题所需的任何材料。这只是一个跟踪溶液和溶质的问题。按照以下步骤完成这个练习：

1. 让我们创建一个函数，以找出在t分钟后的盐含量，给定我们的初始条件：

def salt_content(t):
    salt = 18 #pounds
    brine = 82 #gallons

1. 然后，每分钟都会添加 5 加仑的卤水，其中含有 15 磅（每加仑 3 磅盐每加仑）的盐：

    for i in range(t):
        brine += 5
        salt += 15

1. 现在，每分钟流出 2 加仑的卤水，但其中含有多少盐呢？这要求我们找出每加仑卤水的浓度：

图 12.20：计算每加仑卤水浓度的公式

这可以很容易地转换为代码，如下所示：

        concentration = salt/brine

1. 因此，每分钟离开罐子的盐将是流出的溶液加仑数乘以盐的浓度：

        salt_out = 2*concentration
        salt -= salt_out
        brine -= 2

1. 循环结束后，我们可以打印出卤水和盐的最终数量：

    print(i,brine,salt)

1. 为了解决我们的问题，我们只需运行我们的salt_content函数，t=39：

salt_content(39)

输出如下：

38 199 469.2592152141211

这意味着在 39 分钟后，我们最终得到 469 磅盐。这个数字非常接近解析解，但并不完全相同。我们该怎么做才能得到更准确的结果呢？记住，自然对数的底数e的背后思想是模拟值的恒定变化，而我们只是每分钟计算一次我们溶液的变化。

1. 让我们引入一个名为frac的变量，它将让我们计算分钟的变化：

def salt_content(t,frac=0.001):
    salt = 18 #pounds
    brine = 82 #gallons

1. 参数中的frac=0.001值表示我们将每分钟计算一千次变化。这意味着我们将循环的次数乘以 1,000，或者 1/frac，我们将我们的数量变化乘以frac：

    for i in range(int(t/frac)):
        brine += 5*frac
        salt += 15*frac
        concentration = salt/brine
        salt_out = 2*concentration*frac
        salt -= salt_out
        brine -= 2*frac
    print(i,brine,salt)
salt_content(39)

输出变成了以下内容：

38999 198.99999999966812 470.74539697793307

470.7 磅盐甚至更接近解析解，使用更小的分钟分数并不会改变输出太多。

注意

要访问此特定部分的源代码，请参阅packt.live/2BlX2Tn。

您也可以在packt.live/3dSrEcm上在线运行此示例。

让我们在其他问题上使用这个函数。

练习 12.13：解决混合问题-第 2 部分

一个罐子里含有 10,000 升浓度为每 100 升 1 千克盐的卤水溶液。每秒以 20 升的速度流入罐中含有每 100 升 2 千克盐的卤水。混合物（均匀）以每秒 10 升的速度流出。找出在 5 分钟内罐中有多少盐。按照以下步骤完成这个练习：

1. 因此，我们需要进行一些简单的算术运算来找出我们的初始盐量，但是每 100 升 1 千克盐是 10,000 升中的 100 千克盐，而流入罐中的 20 升中是 0.4 千克盐。这是我们的新函数：

def salt_content(t,frac=.001):
    salt = 100
    brine = 10000
    for i in range(int(t/frac)):
        brine += 20*frac
        salt += 0.4*frac
        concentration = salt/brine
        salt_out = 10*concentration*frac
        salt -= salt_out
        brine -= 10*frac
    return salt

现在，让我们调用salt_content函数：

print(salt_content(5*60))

当我们调用函数时，输出如下：

183.0769053279811

（记住，我们所有的数字都是以秒为单位的，我们想要 5 分钟，因此是5*60参数。）

输出告诉我们，在 5 分钟内溶液中有 183 千克盐。这非常接近解析解。

1. 我们可以通过将硬编码的数字更改为变量来简化我们的任务，因此当我们遇到不同初始卤水量的问题时，例如，我们只需在函数调用中输入不同的数字。我们需要变量来表示初始卤水量（或任何溶液）、溶质的初始量（到目前为止，我们一直在使用盐）、卤水的流入速度、盐的流入速度和卤水的流出速度。以下是如何更改函数的方法：

def salt_content(t,salt_0,brine_0,salt_in,brine_in,v_out,frac=.001):
    salt = salt_0 #pounds
    brine = brine_0 #gallons
    for i in range(int(t/frac)):
        brine += brine_in * frac
        salt += salt_in* frac
        concentration = salt/brine
        salt_out = v_out*concentration* frac
        salt -= salt_out
        brine -= v_out* frac
    return salt

1. 现在，要解决最后一个问题，我们的函数调用将有更多的参数：

salt_content(300,100,10000,0.4,20,10)

输出如下：

183.0769053279811

如您所见，输出应与步骤 1中的相同。让我们将其应用到更多问题上。

注意

要访问此特定部分的源代码，请参阅packt.live/3gkTWOd.

您还可以在packt.live/3eSWF17.上在线运行此示例

练习 12.14：解决混合问题-第 3 部分

一个大桶中含有 100 升糖水混合物，含有 900 克糖。每分钟以 5 克糖每升的速度进入大桶的糖水混合物为 2 升。另一个含有每升 10 克糖的混合物以每分钟 1 升的速度流入大桶。大桶保持混合，所得的混合物以每分钟 3 升的速度从大桶中排出。在 1 小时内找出大桶中的糖量。按照以下步骤完成此练习：

1. 这里唯一的诀窍是总溶液进入速度为每分钟 3 升，总溶质进入速度为每分钟 20 克。以下是函数调用：

salt_content(60,900,100,20,3,3)

1. 输出将如下所示：

705.2374486274181

溶质的量为 705 克。

注意

要访问此特定部分的源代码，请参阅packt.live/2YRWNIl.

您还可以在packt.live/2YRWKfD.上在线运行此示例

练习 12.15：解决混合问题-第 4 部分

如果我们添加纯水会怎样？这会让它更难还是更容易？让我们试试这个。

一个罐子中含有 1200 升水和 18 克盐的卤水混合物。淡水以每分钟 15 升的速度进入罐子，并且罐子被搅拌以保持均匀。一根管子以每分钟 10 升的速度排出混合物。15 分钟后罐子中有多少盐？按照以下步骤完成此练习：

1. 我们可以使用我们的salt_content函数，但盐变量将设置为0。这使得以下函数调用：

print(salt_content(15,18,1200,0,15,10))

1. 15 分钟后的盐含量输出如下：

15.944648402124784

盐含量从 18 克减少到 15.9 克。

注意

访问此特定部分的源代码，请参阅packt.live/2ZsLTIs.

您还可以在packt.live/2AnLrT8.上在线运行此示例

因此，我们已经看到了通常需要大量代数操作才能找到情况的方程的微分方程的几个主题，以便（大概）我们可以插入一个变量并获得所寻找的温度、位置或数量。使用 Python 进行建模和运行模拟已经为我们节省了大量的代数，并且仍然为我们提供了非常准确的答案。

欧拉方法

在大学数学课程中，您学习了所有这些代数方法来求导数和积分以及解决微分方程。我们没有提到拉普拉斯变换，这是解决微分方程的更复杂的方法。现在，关于微分方程的肮脏秘密是，除非您主修工程学，否则学校不会告诉您的是，您在现实生活中遇到的大多数微分方程都没有解析解。

好消息是，数百年来一直有避免混乱代数的数值方法，随着计算机的发明，这些方法已经成为标准。即使存在解析解，数值方法对于实际目的几乎与解析方法一样准确，并且只需花费一小部分时间即可获得解决方案。

欧拉方法的思想非常简单：

1. 从已知点开始。

2. 使用微分方程在此点计算导数。这是曲线在此点处的方向。

3. 向计算出的方向迈出一小步。

4. 重复直到达到所需范围的末尾。

练习 12.16：使用欧拉方法解决微分方程

给定微分方程。您想知道在特定值x处函数y=f(x)的输出。您在图上给出了一个点：f(0) = 1。这意味着，“在每个点上，这个函数的导数是该点的 y 值。”请记住，导数是图上的点朝向或方向。欧拉方法是从初始值开始，即在这种情况下，在(0,1)，并使用微分方程计算到下一个点的方向。微分方程（DE）规定斜率是y值，因此我们在正x方向上迈出一小步：

图 12.21：朝着正确方向迈出小步（希望如此）

导数如下：

图 12.22：函数的导数

因此，ΔY变为以下内容：

图 12.23：计算ΔY 的公式

这是导数和步长的乘积。要找到下一个y值，我们将ΔY添加到先前的y值。在新点上，我们重复这个过程：计算这一点的函数斜率，乘以步长，然后加到当前的y值上。按照以下步骤进行：

1. 让我们编写一个 Python 函数来做到这一点：

def euler(x0,y0,target_x,stepsize):
    x,y = x0,y0
    while x<target_x:
        slope = y #from diff eq
        x += stepsize
        y += stepsize*slope
        print(x,y)
    return y

1. 因此，我们知道初始的x和y。我们想知道x=2时的y；步长可以是½：

print(euler(0,1,2,0.5))

以下是输出：

0.5 1.5
1.0 2.25
1.5 3.375
2.0 5.0625
5.0625

1. 我们不再需要euler函数内的print语句，因此将其注释掉：

        #print(x,y)

1. 第一行是计算斜率的结果，即y值 1，乘以步长½，然后向上移动该距离。如果导数为负，我们将向下移动。在第二行，我们将y值 1.5 乘以步长 0.5，得到 0.75。我们从 0.75 上升到 2.25 等等。在 x 方向上采取小步骤，直到达到目标 x 值 2，我们最终得到y值 5.0625。我们不再需要打印出每一步，但让我们将步长减半 10 次：

for n in [0.5**i for i in range(10)]:
    print(n,euler(0,1,2,n))

以下是输出：

1.0 4.0
0.5 5.0625
0.25 5.9604644775390625
0.125 6.583250172027423
0.0625 6.958666757218805
0.03125 7.166276152788222
0.015625 7.275669793128417
0.0078125 7.3318505987410365
0.00390625 7.3603235532692795
0.001953125 7.374657160341845

因此，步长越小，我们似乎越接近 7.37。这是近似路径的图形：

图 12.24：使用较小的步长获得更好的近似

第四条曲线（右侧的曲线）是我们近似路径的步长为 1 的路径。第三个图形的步长为½，第二条曲线为¼，第一条曲线为 1/8。我们选择微分方程，因为我们知道代数解。

当x为 2 时，e2 = 7.389。添加y=ex*的实际曲线（左侧的第一条曲线），我们可以看到步长越小，近似值越接近实际曲线：

图 12.25：实际曲线添加到第一条曲线的左侧

但是最后的近似值，步长为 0.001953125，需要在 0 和 2 之间进行 1,024 步。很容易理解为什么在计算机发明之前，欧拉方法不如代数方法受欢迎。

注意

要访问此特定部分的源代码，请参阅packt.live/2VEQiaa。

您还可以在packt.live/2ByZvtv上在线运行此示例。

练习 12.17：使用欧拉方法评估函数

在初始值问题（IVP）上使用欧拉方法和步长 0.001：

图 12.26：欧拉方法在初始 VP 上

在这里，y(0) = 1，以便计算近似解y(x)，当x=0.3时：

1. 在euler函数中，在slope=行中输入微分方程：

def euler(x0,y0,target_x,stepsize):
    x,y = x0,y0
    while x<target_x:
        slope = x+y**2 #from diff eq
        x += stepsize
        y += stepsize*slope
    return y

1. 在函数调用中输入适当的参数：

print(euler(0,1,0.3,0.001))

输出应该如下所示：

1.48695561935322

这意味着通过从我们已知的点(0,1)开始，按微分方程指定的方向迈出微小步骤，我们能够预测 1.49 是对应于 x 值 0.3 的近似y值。

注意

要访问此特定部分的源代码，请参阅packt.live/3inHj6S.

您也可以在packt.live/2VFLEbF.上在线运行此示例

Runge-Kutta 方法

由于 Euler 方法仅基于每个点的导数，它存在一个问题，即始终超出或低于真实曲线。毫不奇怪，在 Euler 方法被发明的几个世纪以来，已经对其进行了改进以抵消其缺点。其中一种改进是Runge-Kutta（RK）方法，它将四个近似值平均在一起，其中之一是 Euler 方法，使用区间的开始，另一个使用区间的结束，另外两个近似值使用区间的中点。当这些近似值平均在一起时，中点的近似值被赋予更高的权重。

以下是当 DE 给出时的方程，f(x,y)，起始x和y，x0 和y0，以及步长h：

图 12.27：给出 f(x,y)时的方程

对于下一个y，我们将前面四个近似值平均在一起，k2 和k3 的权重加倍：

图 12.28：对前 4 个近似值进行平均的公式

然后，当然，x增加了h：

图 12.29：将 x 增加 h

这是一大堆代码，但它的功能令人印象深刻。

练习 12.18：实现 Runge-Kutta 方法

在 IVP 上使用 Runge-Kutta 方法和步长 0.2：

图 12.30：步长为 0.2 的 Runge-Kutta 方法

1. 首先，我们定义微分方程。让我们称之为deriv(x,y)：

def deriv(x,y):
    return x**2 + y**2

1. 现在，我们将定义 Runge-Kutta 方法，称之为rk4：

def rk4(x0,y0,target_x,h):
    while x0 <= target_x:
        print(x0,y0)
        k1 = h*deriv(x0,y0)
        k2 = h*deriv(x0 + h/2, y0 + k1/2)
        k3 = h*deriv(x0 + h/2, y0 + k2/2)
        k4 = h*deriv(x0 + h, y0 + k3)
        #These are the values that are fed back into the function:
        y0 = y0 + (1/6)*(k1 + 2*k2 + 2*k3 + k4)
        x0 = x0 + h

1. 当我们从y(0) = 0开始，并且我们想要使用步长为 0.2 来计算y(1)时，这就是我们所说的：

rk4(0,0,1,0.2)

我们的进展如下打印出来：

0 0
0.2 0.0026668666933346665
0.4 0.021360090381533078
0.6 0.0724512003541295
0.8 0.17409018097333867
1.0 0.35025754914481283

1. 使用相同的步长解决相同的问题，但使用 Euler 方法的准确性较低。在euler函数中，将slope=行更改为匹配新微分方程：

        slope = x**2 + y**2

1. 现在，我们使用 Euler 方法打印出解决方案：

print(euler(0,0,1,0.2))

以下是输出：

0.2428567456277198

这与 Runge-Kutta 解决方案并不十分接近。然而，在计算机出现之前，Runge-Kutta 改进可能更有用，因为我们可以简单地减小 Euler 方法中的步长并获得更好的近似值。这是步长为 0.001 的 Euler 方法的相同输出：

print(euler(0,0,1,0.001))

以下是输出：

0.34960542576393877

这只是对用于解方程的数值方法的简要介绍，不是通过代数来解决，而是通过将起始点输入计算机程序并按微分方程指示的方向迈出小步来解决。这是微积分的一个庞大领域，特别是现在免费软件和编程语言，再加上快速的计算机处理器，使以前费力的计算变得轻松。

注意

要访问此特定部分的源代码，请参阅packt.live/3eWxF95.

您也可以在packt.live/3dUlkkg.上在线运行此示例

追踪曲线

微积分中的一个重要话题是追踪曲线，这是一个代理追逐移动目标所经过的路径。由于追逐者直接朝向目标移动，然后目标移动，这种情况导致了各种微分方程。代数可能会变得非常丑陋，这就是微积分教授喜欢这个话题的原因。然而，正如我们所知，微分方程通常是关于寻找一般代数解的，也就是一个函数，而不是一个数字。理论上，我们可以将值代入函数中，以找到特定时间的粒子位置或房间的温度。使用 Python，我们通过对情况进行建模并找到数值解来跳过代数步骤。我们在一般性上失去了一些东西，但在计算的便利性上获得了一些东西。

练习 12.19：找到捕食者捕捉猎物的位置

一只兔子从(0,0)开始，以每秒 1 单位的正y方向奔跑。一只狐狸从(20,0)开始追逐兔子，奔跑速度是兔子的 1.5 倍。狐狸在什么y值处追到兔子？

执行以下步骤：

1. 首先，我们需要从math模块中获取一些函数来测量距离和角度：

from math import sqrt, atan2,sin,cos

1. 我们将编写一个函数，使用毕达哥拉斯定理来测量捕食者位置和猎物位置之间的距离：

def dist(x1,y1,x2,y2):
    """Returns distance from (x1,y1) to (x2,y2)"""
    return sqrt((x1-x2)**2 + (y1-y2)**2)

1. 关键在于猎物和捕食者位置之间y的变化与x的变化代表我们想要的角度的正切。我们知道它们的位置，因此我们使用反正切函数atan2来计算角度，使得捕食者直接指向猎物。我们真正想知道的是如何改变捕食者的x和y坐标，使其朝着猎物移动 1 单位。为了使捕食者朝着猎物转向，我们需要找到两点之间的角度，如下图所示：

图 12.31：捕食者和猎物之间的角度

1. 一旦我们知道变化，我们就可以将向量乘以我们想要的任何速度：

def towards(x1,y1,x2,y2):
    """Returns unit vector in [x,y] format from point
    1 to point 2"""
    dx,dy = x2-x1,y2-y1
    angle = atan2(dy,dx)
    return [cos(angle),sin(angle)]

我们计算x和y的变化，使用arctangent函数计算角度，然后我们使用余弦和正弦来找到捕食者x和y坐标的相应变化，使其朝着猎物走一步。

1. 现在，追逐可以开始了。我们将捕食者和猎物放在它们所在的位置。然后，我们开始一个循环，其中我们将猎物移动一单位（或更准确的增量）：

def chase():
    predator_x,predator_y = 20,0
    predator_v = 1.5 #prey is 1
    prey_x,prey_y = 0,0
    inc = 0.001
    while dist(predator_x,predator_y,prey_x,prey_y) > 0.001:
        prey_y += 1*inc
        p_vec = towards(predator_x,predator_y,\
                        prey_x,prey_y)
        predator_x += predator_v*p_vec[0]*inc
        predator_y += predator_v*p_vec[1]*inc
        #print(dist(predator_x,predator_y,prey_x,prey_y))
    return predator_y

1. 现在，我们运行追逐并打印出捕食者捕捉猎物的y值：

y = chase()
print("Y:",y)
print("dist:",dist(1,1,4,5))
print("towards:",towards(1,1,2,2.732))

输出如下：

Y: 23.997299988652507
dist: 5.0
towards: [0.5000110003630132, 0.8660190526287391]

这非常接近理论值 24。

注意

要访问此特定部分的源代码，请参阅packt.live/3f6x44Z.

您还可以在packt.live/2NO1A7v.上在线运行此示例

练习 12.20：使用乌龟可视化追踪曲线

在这个练习中，我们将可视化捕食者和猎物的路径，这被称为追踪曲线。Python 中有一个内置模块，基于 Logo 编程语言的虚拟乌龟，可以根据我们编写的代码轻松创建可以在屏幕上四处走动的虚拟代理。按照以下步骤完成这个练习：

1. 首先，我们从turtle模块导入函数：

from turtle import *

1. 我们根据所需的左下点设置屏幕的大小，我们将其设置为(-30，-30)，并设置右上点，我们将其设置为(40,40)：

setworldcoordinates(-30,-30,40,40)

1. 设置捕食者和猎物意味着创建一个Turtle对象并设置其颜色、位置和速度。乌龟在行走时会留下路径，所以我们告诉它penup，以防止它在到达起始位置之前绘制。然后，我们告诉它pendown，这样它就会开始绘制：

#set up predator
predator = Turtle()
predator.color("red")
predator.penup()
predator.setpos(20,0)
predator.pendown()
predator.speed(0)

1. 我们通过使乌龟变成绿色并赋予它乌龟的形状来设置猎物：

#set up prey
prey = Turtle()
prey.color("green")
prey.shape("turtle")
prey.setheading(90)
prey.speed(0)

1. pursue函数应该看起来很熟悉，但它有内置函数来计算距离，甚至指向另一个乌龟：

def pursue():
    inc = 0.05
    while predator.distance(prey)>0.05:
        predator.setheading(predator.towards(prey))
        prey.forward(inc)
        predator.forward (1.5*inc)
    print("y:",predator.ycor())

1. 我们将执行pursue函数，然后一旦它打印输出，我们会告诉程序完成，这样图形窗口就不会冻结：

pursue()
done()

1. 如果你运行这个，你可以观看追逐。这是最终输出应该看起来像的：

图 12.32：捕食者的路径是一个对数曲线

1. 扩展：将猎物的路径改为圆形。在使猎物向前移动的行后，添加这行：

prey.left(.3)

这将使猎物每一步左转一小部分度数。但是如果每次转弯都是一样的，最终会形成一个圆。结果路径看起来像这样：

图 12.33：当猎物沿着圆形路径逃离时的追逐曲线

注意

要访问此特定部分的源代码，请参阅packt.live/3dWHDG6。

这一部分目前没有在线交互式示例，需要在本地运行。

位置、速度和加速度

微分方程经常用于研究抛射体的路径，这可以说是微积分的起源。牛顿发明了微积分工具来解决由他对行星运动的研究得出的微分方程，并且表明地球上的自由落体物体受到与轨道行星相同的物理定律的约束。

练习 12.21：计算抛射体离地面的高度

一个球以初速度 29 米/秒向上抛出。它在击中地面之前要多久？按照以下步骤完成这个练习：

1. 在代数课上，我们被引导使用方程来计算抛射体的高度：

图 12.34：计算抛射体高度的公式

这里，h0 是初始高度，v0 是初始向上速度，t是经过的秒数，g是重力加速度，大约 32 英尺或 9.8 米每秒。但是抛射体不使用方程来计算它们的位置；它们只是沿着它们的导数指示的方向运动。

1. 让我们模拟一下：

v = 29
g = 9.8
h = 0
t = 0

因此，对于第一秒，球将以 29 米/秒的速度向上抛出，但会受到每秒 9.8 米的重力减速，这意味着一秒后，它只有29 - 9.8 = 19.2米/秒。因此，一秒后，球应该在空中 19.2 米高。我们每秒重复一次，直到它的高度为 0。

1. height函数应该是这样的：

def height(v0,h0,t):
    """Calculates the height a projectile given the
    initial height and velocity and the elapsed time."""
    v,h = v0,h0
    for i in range(1,t+1):
        v -= g
        h += v
    return h

1. 速度和高度被分配它们的起始值，v0 和h0，然后速度通过g和加速度（由于重力）进行更新，然后高度h通过速度进行更新。我们每秒重复计算并检查球的高度何时返回到零：

for j in range(1,10):
    print(j,round(height(v,h,j),1))

以下是输出：

 –
1 19.2
2 28.6
3 28.2
4 18.0
5 -2.0
6 -31.8
7 -71.4
8 -120.8
9 -180.0

看起来球在 4 到 5 秒之间击中地面。但是当我们将t = 5放入前面的公式中时，我们得到以下结果：

图 12.35：替换计算抛射体高度公式中的值

1. 5 秒后，球应该仍然在空中 22.5 米。我们的代码有什么问题？正如您现在应该知道的那样，球不仅每秒改变一次速度。它的速度不断变化。就像复利一样，我们需要每秒计算新的速度。对于 Python 来说很容易。我们只需引入一个inc变量来增加时间。请记住，这会增加我们循环计算的次数，因此for i in range行也会改变。然后，g和v乘以增量。我们将每半秒重新计算一次：

def height(v0,h0,t):
    """Calculates the height a projectile given the
    initial height and velocity and the elapsed time."""
    inc = 0.5
    v,h = v0,h0
    for i in range(int(t/inc)):
        v -= g*inc
        h += v*inc
    return h

1. 使用相同的代码运行此代码：

for j in range(1,7):
    print(j,round(height(v,h,j),1))

输出现在如下：

1 21.7
2 33.5
3 35.6
4 27.8
5 10.3
6 -17.1

1. 球在空中停留的时间更长，在 5 秒时，它在空中的高度为 10.3 米。如果我们使增量非常小，它应该在 5 秒时更接近 22.5 米。将inc更改为 0.001，再次运行，您将获得以下输出：

1 24.1
2 38.4
3 42.9
4 37.6
5 22.5
6 -2.4

1. 要回答球何时触地的问题，我们将不得不在 5 秒和 6 秒之间进行二分搜索。与以前的搜索一样，我们使用我们的bin_search函数，并更改guess =行以反映我们猜测的数字：

def bin_search(f,lower,upper,target):
    def average(a,b):
        return (a+b)/2
    for i in range(40):
        avg = average(lower,upper)
        guess = f(29,0,avg)
        if guess == target:
            return avg
        if guess < target:
            upper = avg
        else:
            lower = avg
    return avg

1. 我们只需要更改height函数的参数的guess =行。最后一个参数t是我们正在搜索的内容，因此我们正在平均。二分搜索函数将在 5 和 6 之间插入值，并返回返回 0 的t值：

print(bin_search(height,5,6,0))

以下是输出：

5.918000000000575

现在，我们解决二次方程：

29t – 4.9t2 = 0

对于t，我们得到t = 0和t = 5.9184。当然，我们扔出球之前，球的高度为 0，第二个值非常接近我们得到的值。函数的图形如下所示。忽略空气阻力，粒子高度随时间的图形遵循抛物线路径：

图 12.36：没有空气阻力的抛射物路径

这是对我们代码的测试，因为我们有一个很好的公式来检查我们的输出。现在，我们将继续进行更难的关于速度和加速度的微积分问题，其中没有公式可以帮助我们检查答案。

注意

要访问此特定部分的源代码，请参阅packt.live/2VEAkN4。

您还可以在packt.live/2Bzpz7Z上在线运行此示例。

计算带空气阻力的抛射物高度的示例

数学学生被迫研究从代数到微积分的完美抛物线路径中行进的粒子。不幸的是，这不是真实粒子的行进方式。在现实生活中，物体穿过空气或水等介质，并根据介质的密度、其横截面积和其他因素而减速。这导致了一个复杂的施加在抛射物上的力的方程。简而言之，抛射物上的力是由于重力加速度和与其速度的平方成比例的减速。方程如下：

F = mg - kv2

有空气阻力，我们需要知道抛射物的质量m。由于重力加速度g为 9.8 m/s2。变量k是至少三个不同因素的组合，但值k = 0.27对于这种情况产生了现实结果。

与上一个练习一样，我们计算加速度并使用它来更新速度。然后，我们根据速度更新抛射物的位置。

抛射物上的力由两部分组成：由于重力的通常加速度和阻力分量。让我们编写一个 Python 函数来计算：

def force(v,mass,g,k,inc):
    """Returns the downward force on a
    projectile"""
    gravity = mass*g*inc
    drag = k*(v**2)*inc
    if v > 0:
        return gravity + drag
    return gravity - drag

很多时候，我们的值都乘以inc，增量变量，这样我们可以采取更小的步骤来获得更好的近似值，就像我们以前做的那样。gravity和drag变量直接来自力的方程。请注意，如果速度大于 0，抛射物正在向上运动，因此向下的力是重力和阻力力的总和。否则，抛射物正在向下运动，因此重力的力仍然向下，但阻力正在减速，因此我们使用重力和阻力的差异。

现在，我们将调整我们在上一个练习中的height函数，以计算高度等于 0 所需的时间，并添加对我们的force函数的调用：

def height(v0,h0):
    """Calculates the time it takes a projectile given the 
    initial height and velocity to hit the ground."""
    inc = 0.001
    v,h = v0,h0
    t = 0
    while h >= 0:
        v -= force(v,1,9.8,0,inc) #test with k=0
        h += v*inc
        t += inc
    return round(t,1),round(v,1)

在这个函数中，是v -=这一行在起重要作用。速度将受到向下的力的影响。当我们使用k = 0运行时，我们应该得到与上一个问题中相同的时间和结束速度，没有空气阻力：

print(height(29,0))

输出如下：

(5.9, -29.0)

是的；在上一个练习中，抛射物需要 5.9 秒才能到达地面。当没有空气阻力且结束高度与初始高度相同时，结束速度将与初始速度相同，只是方向相反，因此为-29 米/秒。

现在，让我们使用更现实的值k，0.27，看看粒子到达地面需要多长时间，以及它将以多快的速度运行。你有什么预测？

将height函数中的v -=行更改为以下内容：

        v -= force(v,1,9.8,0.27,inc)

当您运行程序时，输出将如下所示：

(2.2, -5.9)

因此，抛射物只用了 2.2 秒就上升并下降，最终速度为-5.9 米/秒。如果我们将具有和没有空气阻力的抛射物的高度图进行对比，我们会发现在有空气阻力的情况下高度要少得多：

图 12.37：带空气阻力的抛射物的高度（内曲线）

这确实是很大的阻力。尝试使用不同的k值，即阻力常数，以获得不同的结束时间和结束速度。这引出了数学和科学中一个非常有趣的概念，即终端速度，当抛射物上下的力相等时，它不再加速。

练习 12.22：计算终端速度

如果您的抛射物从 3000 米的初始高度开始并从飞机上跳下（向下速度为 0），它会达到什么速度？它会简单地继续加速直到抛射物撞击地面吗？

将质量更改为 80 公斤，这是人类的平均体重，k更改为0.27。按照以下步骤完成此练习：

1. 确保您有来自上一个示例的force函数。

2. 修改您的height函数，使其看起来像这样：

def height(v0,h0): 
    """Calculates the velocity of a projectile given the  
    initial height and velocity and the elapsed time.""" 
    inc = 0.001
    v,h = v0,h0 
    t = 0
    for i in range(500): 

1. 这是一个重要的行，我们在其中告诉force函数质量、k的值等等：

        v -= force(v,80,9.8,0.27,inc)
        h += v*inc
        if i % 50 == 0:
            print("v:",round(v,1))
        t += inc

1. 我们进行了 500 次循环，但只在每 50 次循环时打印出速度。让我们用这行来运行它：

height(0,3000)

这是我们收到的输出：

v: -0.8
v: -34.1
v: -48.6
v: -52.6
v: -53.6
v: -53.8
v: -53.9
v: -53.9
v: -53.9
v: -53.9

速度从 0 开始，变得越来越负，直到停止减少。它在大约 54 米/秒左右稳定下来（为负，因为它向下），这大约是每小时 120 英里，人体在自由落体中的终端速度。这是随时间变化的速度图：

图 12.38：带空气阻力的自由落体物体的速度

注意

要访问此特定部分的源代码，请参阅packt.live/2NNmWBM。

您还可以在packt.live/2BUuXCp上在线运行此示例。

现在，让我们完成一个活动，测试我们在本章学到的知识。

活动 12.01：找到粒子的速度和位置

在x-y平面上移动的粒子的速度矢量具有以下分量：

图 12.39：粒子速度矢量的微分方程

找到曲线的切线水平的所有时间（和坐标），然后找到t=1时粒子的速度。

执行以下步骤完成此活动：

1. 编写dx/dt和dy/dt的函数。

2. 循环遍历输出，找到导数为 0 的位置，找到导数从正变负或反之的值。然后，使用二分搜索找到更精确的近似值。

3. 创建一个position函数，并使用循环按照之前给出的导数（位置的变化）来改变粒子的位置，以便在所需的经过时间停止，并打印出x-y坐标。

4. 将您在步骤 2中找到的时间插入position函数中，以找到导数为 0 时粒子的x-y坐标。

5. 您被要求在时间t=1时粒子的速度。使用您得到的微分方程找到粒子速度的垂直和水平分量，并找到以这些分量作为腿的直角三角形的斜边。

注意

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

总结

微积分是一套非常强大的工具，用于模拟真实情况，从热量传递到行星运动。它使我们能够计算函数在瞬间的变化率和复杂曲线下的面积（这些任务仅使用代数和几何的工具似乎是不可能的）。在本章中，我们能够处理值的变化率（导数）作为一个值本身，并使用 Python 循环和函数计算出一些非常精确的结果。模拟导致微分方程的情况，比如抛射物的路径，正是推动第一台电子计算机的发展。

数学课程可能仍然强调对方程的代数解，甚至微分方程，但正如我们在本章中所看到的，使用计算机是模拟现实生活情况的一种简单方法，比如捕食者追逐猎物。我们改变了变量，比如投资中的资金量，混合物中的盐量，以及捕食者的朝向，重复计算了数千次，每一步重新计算数量和距离，得到了非常精确的结果。Python 是设置一些起始条件并让程序运行直到抛射物击中地面或达到终端速度的完美工具。Python 还帮助我们避免了繁琐的代数运算，并让我们通过创建一个简单的模型来蛮力得到答案，比如一个下落的物体或一个捕食者追逐猎物。这是简单的，因为我们不必重复数千次计算——计算机会做。此外，这些数值方法已经用于没有简单代数解的微分方程，甚至适用于那些有解的方程。希望本章已经证明了使用计算机来模拟和分析复杂的现实情况的力量。

现在，您已经学会了如何利用 Python 的循环、变量、条件、函数和列表来解决统计学、概率论和微积分中的复杂问题。您还学会了如何计时执行代码并绘制输出。您使用了 Python 的最先进的数值包numpy来加速计算并操纵数组，适用于各种应用。您还看到 Python 编程被应用于太阳下的每一个数学主题，现在您将能够将其应用于未来遇到的任何现实生活情况。

JMK95

GEA39