Python 数学应用（一）（3）

Python 数学应用（一）（2）https://developer.aliyun.com/article/1506373

曲面和等高线图

Matplotlib 还可以以各种方式绘制三维数据。显示这种数据的两种常见选择是使用表面图或等高线图（类似于地图上的等高线）。在本示例中，我们将看到一种从三维数据绘制表面和绘制三维数据等高线的方法。

准备就绪

要绘制三维数据，需要将其排列成x、y和z分量的二维数组，其中x和y分量必须与z分量的形状相同。为了演示，我们将绘制对应于f(x, y) = x²y³函数的表面。

如何做…

我们想要在-2≤x≤2 和-1≤y≤1 范围内绘制f(x, y) = x²y³函数。第一项任务是创建一个适当的(x, y)对的网格，以便对该函数进行评估：

1. 首先使用np.linspace在这些范围内生成合理数量的点：

X = np.linspace(-2, 2)
Y = np.linspace(-1, 1)

1. 现在，我们需要创建一个网格来创建我们的z值。为此，我们使用np.meshgrid例程：

x, y = np.meshgrid(X, Y)

1. 现在，我们可以创建要绘制的z值，这些值保存了每个网格点上函数的值：

z = x**2 * y**3

1. 要绘制三维表面，我们需要加载一个 Matplotlib 工具箱mplot3d，它随 Matplotlib 包一起提供。这不会在代码中明确使用，但在幕后，它使三维绘图实用程序可用于 Matplotlib：

from mpl_toolkits import mplot3d

1. 接下来，我们创建一个新的图和一组三维坐标轴用于该图：

fig = plt.figure()
ax = fig.add_subplot(projection="3d")  # declare 3d plot

1. 现在，我们可以在这些坐标轴上调用plot_surface方法来绘制数据：

ax.plot_surface(x, y, z)

1. 为三维图添加轴标签非常重要，因为在显示的图上可能不清楚哪个轴是哪个：

ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_zlabel("$z$")

1. 此时我们还应该设置一个标题：

ax.set_title("Graph of the function $f(x) = x²y³$)

您可以使用plt.show例程在新窗口中显示图（如果您在 Python 中交互使用，而不是在 Jupyter 笔记本或 IPython 控制台上使用），或者使用plt.savefig将图保存到文件中。上述序列的结果如下所示：

图 2.6：使用默认设置使用 Matplotlib 生成的三维表面图

1. 等高线图不需要mplot3d工具包，在pyplot接口中有一个contour例程可以生成等高线图。但是，与通常的（二维）绘图例程不同，contour例程需要与plot_surface方法相同的参数。我们使用以下顺序生成绘图：

fig = plt.figure()  # Force a new figure
plt.contour(x, y, z)
plt.title("Contours of $f(x) = x²y³$")
plt.xlabel("$x$")
plt.ylabel("$y$")

结果显示在以下图中：

图 2.7：使用默认设置使用 Matplotlib 生成的等高线图

它是如何工作的…

mplot3d工具包提供了一个Axes3D对象，这是核心 Matplotlib 包中Axes对象的三维版本。当给定projection="3d"关键字参数时，这将被提供给Figure对象上的axes方法。通过在三维投影中在相邻点之间绘制四边形，可以获得表面图。这与用直线连接相邻点来近似二维曲线的方式相同。

plot_surface方法需要提供z值，这些值作为二维数组编码在(x, y)对的网格上的z值。我们创建了我们感兴趣的x和y值的范围，但是如果我们简单地在这些数组的对应值上评估我们的函数，我们将得到一条线上的z值，而不是整个网格上的值。相反，我们使用meshgrid例程，它接受两个X和Y数组，并从中创建一个网格，其中包含X和Y中所有可能的值的组合。输出是一对二维数组，我们可以在其上评估我们的函数。然后我们可以将这三个二维数组全部提供给plot_surface方法。

还有更多…

在前面的部分中描述的例程contour和plot_contour只适用于高度结构化的数据，其中x、y和z分量被排列成网格。不幸的是，现实生活中的数据很少有这么结构化的。在这种情况下，您需要在已知点之间执行某种插值，以近似均匀网格上的值，然后可以绘制出来。执行这种插值的常见方法是通过对(x, y)对的集合进行三角剖分，然后使用每个三角形顶点上的函数值来估计网格点上的值。幸运的是，Matplotlib 有一个方法可以执行所有这些步骤，然后绘制结果，这就是plot_trisurf例程。我们在这里简要解释一下如何使用它：

1. 为了说明plot_trisurf的用法，我们将从以下数据绘制表面和等高线：

x = np.array([ 0.19, -0.82, 0.8 , 0.95, 0.46, 0.71, 
     -0.86, -0.55,   0.75,-0.98, 0.55, -0.17, -0.89,
         -0.4 , 0.48, -0.09, 1., -0.03, -0.87, -0.43])
y = np.array([-0.25, -0.71, -0.88, 0.55, -0.88, 0.23, 
      0.18,-0.06, 0.95, 0.04, -0.59, -0.21, 0.14, 0.94, 
          0.51, 0.47, 0.79, 0.33, -0.85, 0.19])
z = np.array([-0.04, 0.44, -0.53, 0.4, -0.31, 0.13, 
      -0.12, 0.03, 0.53, -0.03, -0.25, 0.03, -0.1 , 
          -0.29, 0.19, -0.03, 0.58, -0.01, 0.55, -0.06])

1. 这次，我们将在同一图中绘制表面和等高线（近似），作为两个单独的子图。为此，我们向包含表面的子图提供projection="3d"关键字参数。我们在三维坐标轴上使用plot_trisurf方法绘制近似表面，并在二维坐标轴上使用tricontour方法绘制近似等高线：

fig = plt.figure(tight_layout=True)  # force new figure
ax1 = fig.add_subplot(1, 2, 1, projection="3d")  # 3d axes
ax1.plot_trisurf(x, y, z)
ax1.set_xlabel("$x$")
ax1.set_ylabel("$y$")
ax1.set_zlabel("$z$")
ax1.set_title("Approximate surface")

1. 现在，我们可以使用以下命令绘制三角剖分表面的等高线：

ax2 = fig.add_subplot(1, 2, 2)  # 2d axes
ax2.tricontour(x, y, z)
ax2.set_xlabel("$x$")
ax2.set_ylabel("$y$")
ax2.set_title("Approximate contours")

我们在图中包含tight_layout=True关键字参数，以避免稍后调用plt.tight_layout例程。结果如下所示：

图 2.8：使用三角剖分生成的近似表面和等高线图

除了表面绘图例程外，Axes3D对象还有一个用于简单三维绘图的plot（或plot3D）例程，其工作方式与通常的plot例程完全相同，但在三维坐标轴上。该方法还可用于在其中一个轴上绘制二维数据。

自定义三维图

等高线图可能会隐藏表示的表面的一些细节，因为它们只显示“高度”相似的地方，而不显示值是多少，甚至与周围的值相比如何。在地图上，这可以通过在特定等高线上打印高度来解决。表面图更具启发性，但是将三维对象投影到二维以在屏幕上显示可能会模糊一些细节。为了解决这些问题，我们可以自定义三维图（或等高线图）的外观，以增强图表并确保我们希望突出显示的细节清晰可见。最简单的方法是通过更改图表的颜色映射来实现这一点。

在这个示例中，我们将使用binary颜色映射的反转。

准备工作

我们将为以下函数生成表面图：

我们生成应该绘制的点，就像在前一个示例中一样：

X = np.linspace(-2, 2)
Y = np.linspace(-2, 2)
x, y = np.meshgrid(X, Y)
t = x**2 + y**2  # small efficiency
z = np.cos(2*np.pi*t)*np.exp(-t)

如何做…

Matplotlib 有许多内置的颜色映射，可以应用于图表。默认情况下，表面图是用根据光源进行着色的单一颜色绘制的（请参阅本示例的更多信息部分）。颜色映射可以显著改善图表的效果。以下步骤显示了如何向表面和等高线图添加颜色映射：

1. 首先，我们只需应用内置的颜色映射之一binary_r，通过向plot_surface例程提供cmap="binary_r"关键字参数来实现：

fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.plot_surface(x, y, z, cmap="binary_r")
ax.set_title("Surface with colormap")
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_zlabel("$z$")

结果是一个图（图 2.9），其中表面的颜色根据其值而变化，颜色映射的两端具有最极端的值——在本例中，z值越大，灰度越浅。请注意，下图中的不规则性是由网格中相对较少的点造成的：

图 2.9：应用灰度颜色映射的表面图

颜色映射适用于表面绘图以外的其他绘图类型。特别是，颜色映射可以应用于等高线图，这有助于区分代表较高值和代表较低值的等高线。

1. 对于等高线图，更改颜色映射的方法是相同的；我们只需为cmap关键字参数指定一个值：

fig = plt.figure()
plt.contour(x, y, z, cmap="binary_r")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.title("Contour plot with colormap set")

上述代码的结果如下所示：

图 2.10：具有替代颜色映射设置的等高线图

图中较深的灰色阴影对应于 z 的较低值。

工作原理…

颜色映射通过根据比例尺分配 RGB 值来工作——颜色映射。首先，对值进行归一化，使其介于0和1之间，通常通过线性变换来实现，将最小值取为0，最大值取为1。然后将适当的颜色应用于表面绘图的每个面（或者在另一种类型的绘图中是线）。

Matplotlib 带有许多内置的颜色映射，可以通过简单地将名称传递给cmap关键字参数来应用。这些颜色映射的列表在文档中给出（matplotlib.org/tutorials/colors/colormaps.html），还有一个反转的变体，通过在所选颜色映射的名称后添加_r后缀来获得。

还有更多…

应用颜色映射中的归一化步骤是由Normalize类派生的对象执行的。Matplotlib 提供了许多标准的归一化例程，包括LogNorm和PowerNorm。当然，您也可以创建自己的Normalize子类来执行归一化。可以使用plot_surface或其他绘图函数的norm关键字添加替代Normalize子类。

对于更高级的用途，Matplotlib 提供了一个接口，用于使用光源创建自定义阴影。这是通过从matplotlib.colors包中导入LightSource类，然后使用该类的实例根据z值对表面元素进行着色来完成的。这是使用LightSource对象的shade方法完成的：

from matplotlib.colors import LightSource
light_source = LightSource(0, 45)  # angles of lightsource
cmap = plt.get_cmap("binary_r")
vals = light_source.shade(z, cmap)
surf = ax.plot_surface(x, y, z, facecolors=vals)

如果您希望了解更多关于这个工作原理的内容，可以在 Matplotlib 库中查看完整的示例。

进一步阅读

Matplotlib 包非常庞大，我们在这么短的篇幅内几乎无法充分展现它。文档中包含的细节远远超过了这里提供的内容。此外，还有一个大型的示例库（matplotlib.org/gallery/index.html#），其中包含了比本书中更多的包功能。

还有其他构建在 Matplotlib 之上的包，为特定应用程序提供了高级绘图方法。例如，Seaborn 库提供了用于可视化数据的例程（seaborn.pydata.org/）。

第四章：微积分和微分方程

在本章中，我们将讨论与微积分相关的各种主题。微积分是涉及微分和积分过程的数学分支。从几何上讲，函数的导数代表函数曲线的梯度，函数的积分代表函数曲线下方的面积。当然，这些特征只在某些情况下成立，但它们为本章提供了一个合理的基础。

我们首先来看一下简单函数类的微积分：多项式。在第一个示例中，我们创建一个表示多项式的类，并定义不同 iate 和积分多项式的方法。多项式很方便，因为多项式的导数或积分再次是多项式。然后，我们使用 SymPy 包对更一般的函数进行符号微分和积分。之后，我们看到使用 SciPy 包解方程的方法。接下来，我们将注意力转向数值积分（求积）和解微分方程。我们使用 SciPy 包来解常微分方程和常微分方程组，然后使用有限差分方案来解简单的偏微分方程。最后，我们使用快速傅里叶变换来处理嘈杂的信号并滤除噪音。

在本章中，我们将涵盖以下示例：

• 使用多项式和微积分
  
• 使用 SymPy 进行符号微分和积分
  
• 解方程
  
• 使用 SciPy 进行数值积分
  
• 使用数值方法解简单的微分方程
  
• 解微分方程组
  
• 使用数值方法解偏微分方程
  
• 使用离散傅里叶变换进行信号处理
  

技术要求

除了科学 Python 包 NumPy 和 SciPy 之外，我们还需要 SymPy 包。可以使用您喜欢的软件包管理器（如pip）来安装它：

python3.8 -m pip install sympy

本章的代码可以在 GitHub 存储库的Chapter 03文件夹中找到，网址为github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2003。

查看以下视频以查看代码的实际操作：bit.ly/32HuH4X。

使用多项式和微积分

多项式是数学中最简单的函数之一，定义为一个求和：

x代表要替换的占位符，a[i]是一个数字。由于多项式很简单，它们提供了一个很好的方式来简要介绍微积分。微积分涉及函数的微分和积分。积分大致上是反微分，因为先积分然后微分会得到原始函数。

在本示例中，我们将定义一个表示多项式的简单类，并为该类编写方法以执行微分和积分。

准备工作

从几何上讲，通过微分得到的导数是函数的梯度，通过积分得到的积分是函数曲线与x轴之间的面积，考虑到曲线是在轴的上方还是下方。在实践中，微分和积分是通过一组规则和特别适用于多项式的标准结果来进行符号化处理。

本示例不需要额外的软件包。

如何做…

以下步骤描述了如何创建表示多项式的类，并为该类实现微分和积分方法：

1. 让我们首先定义一个简单的类来表示多项式：

class Polynomial:
    """Basic polynomial class"""
    def __init__(self, coeffs):
        self.coeffs = coeffs
    def __repr__(self):
        return f"Polynomial({repr(self.coeffs)})"
    def __call__(self, x):
        return sum(coeff*x**i for i, coeff 
                in enumerate(self.coeffs))

1. 现在我们已经为多项式定义了一个基本类，我们可以继续实现这个Polynomial类的微分和积分操作，以说明这些操作如何改变多项式。我们从微分开始。我们通过将当前系数列表中的每个元素（不包括第一个元素）相乘来生成新的系数。我们使用这个新的系数列表来创建一个新的Polynomial实例，然后返回：

def differentiate(self):
        """Differentiate the polynomial and return the derivative"""
        coeffs = [i*c for i, c in enumerate(self.coeffs[1:], start=1)]
        return Polynomial(coeffs)

1. 要实现积分方法，我们需要创建一个包含由参数给出的新常数（转换为浮点数以保持一致性）的新系数列表。然后我们将旧系数除以它们在列表中的新位置，加到这个系数列表中：

def integrate(self, constant=0):
        """Integrate the polynomial, returning the integral"""
        coeffs = [float(constant)] 
        coeffs += [c/i for i, c in enumerate(self.coeffs, start=1)]
        return Polynomial(coeffs)

1. 最后，为了确保这些方法按预期工作，我们应该用一个简单的案例测试这两种方法。我们可以使用一个非常简单的多项式来检查，比如x² - 2x + 1：

p = Polynomial([1, -2, 1])
p.differentiate()
# Polynomial([-2, 2])
p.integrate(constant=1)
# Polynomial([1.0, 1.0, -1.0, 0.3333333333])

工作原理…

多项式为我们提供了一个对微积分基本操作的简单介绍，但对于其他一般类的函数来说，构建 Python 类并不那么容易。也就是说，多项式非常有用，因为它们被很好地理解，也许更重要的是，对于多项式的微积分非常容易。对于变量x的幂，微分的规则是乘以幂并减少 1，因此x^n*变为*nx(n-1)。

积分更复杂，因为函数的积分不是唯一的。我们可以给积分加上任意常数并得到第二个积分。对于变量x的幂，积分的规则是将幂增加 1 并除以新的幂，因此x^n*变为*x(n+1)/(n+1)，所以要对多项式进行积分，我们将每个x的幂增加 1，并将相应的系数除以新的幂。

我们在这个食谱中定义的Polynomial类相当简单，但代表了核心思想。多项式由其系数唯一确定，我们可以将其存储为一组数值值的列表。微分和积分是我们可以对这个系数列表执行的操作。我们包括一个简单的__repr__方法来帮助显示Polynomial对象，以及一个__call__方法来促进在特定数值上的评估。这主要是为了演示多项式的评估方式。

多项式对于解决涉及评估计算昂贵函数的某些问题非常有用。对于这样的问题，我们有时可以使用某种多项式插值，其中我们将一个多项式“拟合”到另一个函数，然后利用多项式的性质来帮助解决原始问题。评估多项式比原始函数要“便宜”得多，因此这可能会大大提高速度。这通常是以一些精度为代价的。例如，辛普森法则用二次多项式逼近曲线下的面积，这些多项式是由三个连续网格点定义的间隔内的。每个二次多项式下面的面积可以通过积分轻松计算。

还有更多…

多项式在计算编程中扮演的角色远不止是展示微分和积分的效果。因此，NumPy 包中提供了一个更丰富的Polynomial类，numpy.polynomial。NumPy 的Polynomial类和各种派生子类在各种数值问题中都很有用，并支持算术运算以及其他方法。特别是，有用于将多项式拟合到数据集合的方法。

NumPy 还提供了从Polynomial派生的类，表示各种特殊类型的多项式。例如，Legendre类表示一种特定的多项式系统，称为Legendre多项式。Legendre 多项式是为满足*-1 ≤ x ≤ 1的x*定义的，并且形成一个正交系统，这对于诸如数值积分和有限元方法解决偏微分方程等应用非常重要。Legendre 多项式使用递归关系定义。我们定义

对于每个n ≥ 2，我们定义第n个 Legendre 多项式满足递推关系，

还有一些所谓的正交（系统的）多项式，包括Laguerre多项式*，Chebyshev 多项式和Hermite 多项式。

参见

微积分在数学文本中有很好的文档记录，有许多教科书涵盖了从基本方法到深层理论的内容。正交多项式系统在数值分析文本中也有很好的文档记录。

使用 SymPy 进行符号微分和积分

在某些时候，您可能需要区分一个不是简单多项式的函数，并且可能需要以某种自动化的方式来做这件事，例如，如果您正在编写教育软件。Python 科学堆栈包括一个名为 SymPy 的软件包，它允许我们在 Python 中创建和操作符号数学表达式。特别是，SymPy 可以执行符号函数的微分和积分，就像数学家一样。

在这个示例中，我们将创建一个符号函数，然后使用 SymPy 库对这个函数进行微分和积分。

准备工作

与其他一些科学 Python 软件包不同，文献中似乎没有一个标准的别名来导入 SymPy。相反，文档在几个地方使用了星号导入，这与 PEP8 风格指南不一致。这可能是为了使数学表达更自然。我们将简单地按照其名称sympy导入模块，以避免与scipy软件包的标准缩写sp混淆（这也是sympy的自然选择）：

import sympy

在这个示例中，我们将定义一个表示函数的符号表达式

如何做…

使用 SymPy 软件包进行符号微分和积分（就像您手工操作一样）非常容易。按照以下步骤来看看它是如何完成的：

1. 一旦导入了 SymPy，我们就定义将出现在我们的表达式中的符号。这是一个没有特定值的 Python 对象，就像数学变量一样，但可以在公式和表达式中表示许多不同的值。对于这个示例，我们只需要定义一个符号用于x，因为除此之外我们只需要常数（文字）符号和函数。我们使用sympy中的symbols例程来定义一个新符号。为了保持符号简单，我们将这个新符号命名为x：

x = sympy.symbols('x')

1. 使用symbols函数定义的符号支持所有算术运算，因此我们可以直接使用我们刚刚定义的符号x构造表达式：

f = (x**2 - 2*x)*sympy.exp(3 - x)

1. 现在我们可以使用 SymPy 的符号微积分能力来计算f的导数，即对f进行微分。我们使用sympy中的diff例程来完成这个操作，它对指定的符号进行符号表达式微分，并返回导数的表达式。这通常不是以最简形式表达的，因此我们使用sympy.simplify例程来简化结果：

fp = sympy.simplify(sympy.diff(f))  # (x*(2 - x) + 2*x - 2)
                                                       *exp(3 - x)

1. 我们可以通过以下方式检查使用 SymPy 进行符号微分的结果是否正确，与手工计算的导数相比，定义为 SymPy 表达式：

fp2 = (2*x - 2)*sympy.exp(3 - x) - (x**2 - 2*x)*sympy.exp(3 - x)

1. SymPy 相等性测试两个表达式是否相等，但不测试它们是否在符号上等价。因此，我们必须首先简化我们希望测试的两个语句的差异，并测试是否等于0：

sympy.simplify(fp2 - fp) == 0  # True

1. 我们可以使用 SymPy 通过integrate函数对函数f进行积分。还可以通过将其作为第二个可选参数提供来提供要执行积分的符号：

F = sympy.integrate(f, x)  # -x**2*exp(3 - x)

它是如何工作的…

SymPy 定义了表示某些类型表达式的各种类。例如，由Symbol类表示的符号是原子表达式的例子。表达式是以与 Python 从源代码构建抽象语法树的方式构建起来的。然后可以使用方法和标准算术运算来操作这些表达式对象。

SymPy 还定义了可以在Symbol对象上操作以创建符号表达式的标准数学函数。最重要的特性是能够执行符号微积分 - 而不是我们在本章剩余部分中探索的数值微积分 - 并给出对微积分问题的精确（有时称为解析）解决方案。

SymPy 软件包中的diff例程对这些符号表达式进行微分。这个例程的结果通常不是最简形式，这就是为什么我们在配方中使用简化例程来简化导数的原因。integrate例程用给定的符号对scipy表达式进行符号积分。（diff例程还接受一个符号参数，用于指定微分的符号。）这将返回一个其导数为原始表达式的表达式。这个例程不会添加积分常数，这在手工积分时是一个好的做法。

还有更多…

SymPy 可以做的远不止简单的代数和微积分。它有各种数学领域的子模块，如数论、几何和其他离散数学（如组合数学）。

SymPy 表达式（和函数）可以构建成 Python 函数，可以应用于 NumPy 数组。这是使用sympy.utilities模块中的lambdify例程完成的。这将 SymPy 表达式转换为使用 SymPy 标准函数的 NumPy 等价函数来数值评估表达式。结果类似于定义 Python Lambda，因此得名。例如，我们可以使用这个例程将这个配方中的函数和导数转换为 Python 函数：

from sympy.utilities import lambdify
lam_f = lambdify(x, f)
lam_fp = lambdify(x, fp)

lambdify例程接受两个参数。第一个是要提供的变量，上一个代码块中的x，第二个是在调用此函数时要评估的表达式。例如，我们可以评估之前定义的 lambdified SymPy 表达式，就好像它们是普通的 Python 函数一样：

lam_f(4)  # 2.9430355293715387
lam_fp(7)  # -0.4212596944408861

我们甚至可以在 NumPy 数组上评估这些 lambdified 表达式：

lam_f(np.array([0, 1, 2]))  # array([ 0\. , -7.3890561, 0\. ])

lambdify例程使用 Python 的exec例程来执行代码，因此不应该与未经过消毒的输入一起使用。

解方程

许多数学问题最终归结为解形式为f(x) = 0 的方程，其中f是单变量函数。在这里，我们试图找到一个使方程成立的x的值。使方程成立的x的值有时被称为方程的根。有许多算法可以找到这种形式的方程的解。在这个配方中，我们将使用牛顿-拉弗森和弦截法来解决形式为f(x) = 0 的方程。

牛顿-拉弗森方法（牛顿法）和弦截法是良好的标准根查找算法，几乎可以应用于任何情况。这些是迭代方法，从一个根的近似值开始，并迭代改进这个近似值，直到它在给定的容差范围内。

为了演示这些技术，我们将使用使用 SymPy 进行符号计算配方中定义的函数

它对所有实数x都有定义，并且恰好有两个根，一个在x=0，另一个在x=2。

准备工作

SciPy 包含用于解方程的例程（以及许多其他内容）。根查找例程可以在scipy包的optimize模块中找到。

如果你的方程不是形式上的f(x) = 0，那么你需要重新排列它，使其成为这种情况。这通常不太困难，只需要将右侧的任何项移到左侧即可。例如，如果你希望找到函数的不动点，也就是当g(x)= x时，我们会将方法应用于由f(x) =g(x)*- x.*给出的相关函数。

如何操作…

optimize包提供了用于数值根查找的例程。以下说明描述了如何使用该模块中的newton例程：

1. optimize模块没有列在scipy命名空间中，所以你必须单独导入它：

from scipy import optimize

1. 然后我们必须在 Python 中定义这个函数及其导数：

from math import exp
def f(x):
    return x*(x - 2)*exp(3 - x)

1. 这个函数的导数在前一个配方中被计算出来：

def fp(x):
    return -(x**2 - 4*x + 2)*exp(3 - x)

1. 对于牛顿-拉弗森和割线法，我们使用optimize中的newton例程。割线法和牛顿-拉弗森法都需要函数和第一个参数以及第一个近似值x0作为第二个参数。要使用牛顿-拉弗森法，我们必须提供f的导数，使用fprime关键字参数：

optimize.newton(f, 1, fprime=fp) # Using the Newton-Raphson method
# 2.0

1. 使用割线法时，只需要函数，但是我们必须提供根的前两个近似值；第二个作为x1关键字参数提供：

optimize.newton(f, 1., x1=1.5) # Using x1 = 1.5 and the secant method
# 1.9999999999999862

牛顿-拉弗森法和割线法都不能保证收敛到根。完全有可能方法的迭代只是在一些点之间循环（周期性）或者波动剧烈（混沌）。

工作原理…

对于具有导数f’(x)和初始近似值x[0]的函数f(x)，牛顿-拉弗森方法使用以下公式进行迭代定义

对于每个整数i ≥0。从几何上讲，这个公式是通过考虑梯度的负方向（所以函数是递减的）如果f(x[i])>0或正方向（所以函数是递增的）如果f(x[i]) <o。

割线法基于牛顿-拉弗森法，但是用近似值替换了一阶导数

当x[i]-x[i-1]足够小时，这意味着方法正在收敛，这是一个很好的近似值。不需要函数f的导数的代价是我们需要一个额外的初始猜测来启动方法。该方法的公式如下

一般来说，如果任一方法得到一个足够接近根的初始猜测（割线法的猜测），那么该方法将收敛于该根。牛顿-拉弗森法在迭代中导数为零时也可能失败，此时公式未被很好地定义。

还有更多…

本配方中提到的方法是通用方法，但在某些情况下可能有更快或更准确的方法。广义上讲，根查找算法分为两类：在每次迭代中使用函数梯度信息的算法（牛顿-拉弗森、割线、Halley）和需要根位置的界限的算法（二分法、regula-falsi、Brent）。到目前为止讨论的算法属于第一类，虽然通常相当快，但可能无法收敛。

第二种算法是已知根存在于指定区间内的算法a ≤**x ≤b。我们可以通过检查f(a)和f(b)是否有不同的符号来检查根是否在这样的区间内，也就是说，f(a) <0<f(b)或f(b) <0<f(a)其中一个为真。（当然，前提是函数是连续的，在实践中往往是这样。）这种类型的最基本算法是二分法算法，它重复地将区间二分，直到找到根的足够好的近似值。基本前提是在a和b之间分割区间，并选择函数改变符号的区间。该算法重复，直到区间非常小。以下是 Python 中此算法的基本实现：

from math import copysign
def bisect(f, a, b, tol=1e-5):
    """Bisection method for root finding"""
    fa, fb = f(a), f(b)
    assert not copysign(fa, fb) == fa, "Function must change signs"
    while (b - a) > tol:
        m = (b - a)/2 # mid point of the interval
        fm = f(m)
        if fm == 0:
            return m
        if copysign(fm, fa) == fm: # fa and fm have the same sign
            a = m
            fa = fm
        else: # fb and fm have the same sign
            b = m
        return a

该方法保证收敛，因为在每一步中，距离b-a减半。但是，可能需要比牛顿-拉弗森或割线法更多的迭代次数。optimize中也可以找到二分法的版本。这个版本是用 C 实现的，比这里呈现的版本要高效得多，但是在大多数情况下，二分法不是最快的方法。

Brent 的方法是对二分法的改进，在optimize模块中作为brentq可用。它使用二分和插值的组合来快速找到方程的根：

optimize.brentq(f, 1.0, 3.0)  # 1.9999999999998792

重要的是要注意，涉及括号（二分法、regula-falsi、Brent）的技术不能用于找到复变量的根函数，而不使用括号（Newton、割线、Halley）的技术可以。

使用 SciPy 进行数值积分

积分可以解释为曲线与x轴之间的区域，根据这个区域是在轴的上方还是下方进行标记。有些积分无法直接使用符号方法计算，而必须进行数值近似。其中一个经典例子是高斯误差函数，在第一章的基本数学函数部分中提到。这是由以下公式定义的

这里出现的积分也无法通过符号方法计算。

在本示例中，我们将看到如何使用 SciPy 包中的数值积分例程来计算函数的积分。

准备工作

我们使用scipy.integrate模块，其中包含几个用于计算数值积分的例程。我们将此模块导入如下：

from scipy import integrate

操作步骤…

以下步骤描述了如何使用 SciPy 进行数值积分：

1. 我们评估出现在误差函数定义中的积分在x = 1处的值。为此，我们需要在 Python 中定义被积函数（出现在积分内部）：

def erf_integrand(t):
    return np.exp(-t**2)

scipy.integrate中有两个主要例程用于执行数值积分（求积），可以使用。第一个是quad函数，它使用 QUADPACK 执行积分，第二个是quadrature。

1. quad例程是一个通用的积分工具。它期望三个参数，即要积分的函数（erf_integrand），下限（-1.0）和上限（1.0）：

val_quad, err_quad = integrate.quad(erf_integrand, -1.0, 1.0)
# (1.493648265624854, 1.6582826951881447e-14)

第一个返回值是积分的值，第二个是误差的估计。

1. 使用quadrature例程重复计算，我们得到以下结果。参数与quad例程相同：

val_quadr, err_quadr = integrate.quadrature(erf_integrand, -1.0,
    1.0)
# (1.4936482656450039, 7.459897144457273e-10)

输出与代码的格式相同，先是积分的值，然后是误差的估计。请注意，quadrature例程的误差更大。这是因为一旦估计的误差低于给定的容差，方法就会终止，当调用例程时可以修改这个容差。

工作原理…

大多数数值积分技术都遵循相同的基本过程。首先，我们选择积分区域中的点x[i]，对于i = 1, 2,…, n，然后使用这些值和值f(x[i])来近似积分。例如，使用梯形法则，我们通过以下方式近似积分

其中a < x[1]< x[2]< … < x[n-1]< b，h是相邻x[i]值之间的（公共）差异，包括端点a和b。这可以在 Python 中实现如下：

def trapezium(func, a, b, n_steps):
    """Estimate an integral using the trapezium rule"""
    h = (b - a) / n_steps
    x_vals = np.arange(a + h, b, h) 
    y_vals = func(x_vals)
    return 0.5*h*(func(a) + func(b) + 2.*np.sum(y_vals))

quad和quadrature使用的算法比这复杂得多。使用这个函数来近似使用trapezium积分erf_integrand的积分得到的结果是 1.4936463036001209，这与quad和quadrature例程的近似结果在 5 位小数的情况下是一致的。

quadrature例程使用固定容差的高斯积分，而quad例程使用 Fortran 库 QUADPACK 例程中实现的自适应算法。对两个例程进行计时，我们发现对于配方中描述的问题，quad例程大约比quadrature例程快 5 倍。quad例程在大约 27 微秒内执行，平均执行 1 百万次，而quadrature例程在大约 134 微秒内执行。（您的结果可能会因系统而异。）

还有更多…

本节提到的例程需要知道被积函数，但情况并非总是如此。相反，可能是我们知道一些(x,y)对，其中y = f(x)，但我们不知道要在额外点上评估的函数f。在这种情况下，我们可以使用scipy.integrate中的采样积分技术之一。如果已知点的数量非常大，并且所有点都是等间距的，我们可以使用 Romberg 积分来很好地近似积分。为此，我们使用romb例程。否则，我们可以使用梯形法则的变体（如上所述）使用trapz例程，或者使用simps例程使用辛普森法则。

数值解简单微分方程

微分方程出现在一个数量根据给定关系演变的情况中，通常是随着时间的推移。它们在工程和物理学中非常常见，并且自然地出现。一个经典的（非常简单）微分方程的例子是牛顿提出的冷却定律。物体的温度以与当前温度成比例的速率冷却。从数学上讲，这意味着我们可以写出物体在时间t > 0时的温度T的导数，使用微分方程

常数k是一个确定冷却速率的正常数。这个微分方程可以通过首先“分离变量”，然后积分和重新排列来解析地解决。执行完这个过程后，我们得到了一般解

其中*T[0]*是初始温度。

在这个配方中，我们将使用 SciPy 的solve_ivp例程数值地解决一个简单的常微分方程。

Python 数学应用（一）（4）https://developer.aliyun.com/article/1506376