用Python学《微积分B》（Newton-Leibniz公式）

  微积分的出现源于极限论，而极限的描述借助于“不等式”。“微分”（导数）和“积分”（定积分）一开始只是两个种重要的“极限问题”：微分学研究的是“函数的变化量与自变量的变化量比值的极限”问题，几何上描述为“函数曲线在某点的切线的斜率”问题，其中这个切线就是割线取极限得到的；而积分学研究的是如何“在某个区间内用无限小分割的矩形的面积和来近似函数曲线与x轴围成的面积”，即“Riemann和”，需要注意的是，“Riemann和”是一个极限值。从公式上来看，如下：

df(x)dx=limΔx→0f(x+Δx)−f(x)Δx∫baf(x)dx=limλ→0+∑i=0n[f(ξi)∗Δxi]

$$\frac{\mathrm{d} f(x)}{\mathrm{d} x}=\lim_{\Delta x\rightarrow 0}\frac{f(x+\Delta x)-f(x)}{\Delta x} \\ \int_{a}^{b}f(x)\mathrm{d}x=\lim_{\lambda \rightarrow 0^{+}}\sum_{i=0}^{n}[f(\xi _{i})*\Delta x_{i}]$$
  既然微分和积分都是极限问题，那么它们之间是否有联系呢？答案是肯定的。天才数学家们在微分和积分这两条平行发展的“极限问题研究”的道路之间架起了一座桥梁，那就是“Newton-Leibniz公式”。本文主要分享我在学习《微积分B》7-3节“Newton-Leibniz公式”的体会，并使用pyhon对一些知识点进行扩展和辅助求解课后练习题。
注：除了Newton-Leibniz公式，微分中值定理也是微积分中的一座桥梁，它建立了“函数”与“导数”之间的不严格定量关系（Newton-Leibniz公式是严格定量关系），特别地，在极限情况下非常有用。

一、变上限积分

  要讲Newton-Leibniz公式，必须先讲“变上限积分”，因为Newton-Leibniz公式是在研究“变上限积分”的时候提炼出来的，也是以“变限积分”的形式来描述的。
1，变上限积分不是积分
  变上限积分这个名字有很大的欺骗性，它本质是一个函数，而不是一个积分。它描述的是：对于给定函数（曲线），它在给定左区间（左端点）时的积分值与右端点取值的对应关系（映射）。如下图所示：

  形象点说，就是上图中的面积与x的坐标之间的对应关系（映射），我们可以把这个对应关系记作$S=F(x)$，则有：

F(x)=S=∫xaf(t)dt

$$F(x)=S=\int_a^x f(t)\mathrm{d}t$$
  这就是变上限积分。很明显，它本质上是一个函数，只是这个函数比较特别，它的表达式中含有积分符号，而且自变量是积分上限。
注：上图画图代码参考matplotlib示例
2，定积分的路径和方向
  有了变上限积分，很容易能想到变下限积分，甚至上下限都是x的函数的“变限积分”。其中：
$$\int_x^a f(t)\mathrm{d}t=-\int_a^x f(t)\mathrm{d}t$$
  这个公式看似简单，却不太好证明，也不能直接从“面积法”上感受到。原课程是应用定积分的性质之一$\int_b^a f(x)\mathrm{d}x=-\int_a^b f(x)\mathrm{d}x$来证明的，我觉得不太彻底，因为这条性质本身就让人疑惑。可惜我没有想到更彻底的证明办法，但是我想到了一个直观的现象来辅助理解：
  “Math is fun”上用“水龙头与水缸”（Tap and Tank）的例子来演示曲线与曲线的积分：

  $如上图，\int_a^b f(x)\mathrm{d}x$表示时间正向流动时，我们观察到的水缸中的水随时间的变化情况——水量越来越多。如果时间倒流呢？或者说我们把影片倒过来放映，你会发现水缸中的水越来越少了。而这个时间倒流，就相当于$\int_b^a f(x)\mathrm{d}x$。
  很显然，积分的方向相反（时间流向），积分结果相反，但绝对值相等。更进一步，我们可以想象一下：如果时间变快了，会产生什么现象？是不是和使用快进模式看电影的感觉一样？甚至时间流逝不是恒速而是变速的时候，又会怎么样？
3，上限变量的复合
  有了上面时间流动的例子，再来理解变上限积分的上限是x的函数时的表达式就相对容易了。比如：
$$F(x)=\int_c^{g(x)} f(t)=\int_{g(a)}^{g(x)} f(t)\mathrm{d}t$$
  我还是以“看电影”为例：看电影的时候，同时存在两个时间标尺（坐标轴），电影中的时间尺（t轴）和现实中的时间尺（x轴）。如果我们选择正常播放，那么这两个时间尺（坐标轴）的映射关系是$x(t)=t$，如果我们选择2倍快进播放，那么这个映射关系就是$x(t)=2t$，我们可以发挥想象，建立任意的$x(t)=g(t)$的映射关系。将上面的公式进行拆分：
$$\left\{\begin{matrix} y = F(x)&\cdot \cdot \cdot (1) \\ \; \;y=\int_{g(a)}^{g(x)} f(t)\mathrm{d}t& \cdot \cdot \cdot (2) \end{matrix}\right.$$
  将(1)式和(2)式分别画成图，如下：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import sympy as sp

u = sp.Symbol('u')
b = 16   # interval

# $f=\sqrt{1-(\frac{t}{4})^{2}}$
def func(t):
    return np.sqrt(1 - (t / b) ** 2)

# integral of f
def inte(x):
    return sp.integrate(sp.sqrt(1 - (u / b) ** 2), (u, 0, x ** 2))

a, x = 1, 3.6  # integral limits
t = np.linspace(0, b)
y = func(t)
t2 = np.sqrt(t)
y2 = [inte(n) for n in t2]
# plot curve
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4))
plt.plot(t, y, 'r', linewidth=2)

plt.ylim(ymin=0)
axes[0].plot(t2, y2, 'g', linewidth=2)
axes[0].text(0.5 * (a + x), 5, r"$y=F(x)$",
         horizontalalignment='center', fontsize=20)
axes[1].text(0.5 * (a + x ** 2), 0.4, r"$S=\int_a^{x^2} \sqrt{1-(\frac{t}{4})^{2}}\mathrm{d}t$",
         horizontalalignment='center', fontsize=20)

# Make the shaded region
it = np.linspace(a, x ** 2)
iy = func(it)
verts = [(a, 0)] + list(zip(it, iy)) + [(x ** 2, 0)]
poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')
axes[1].add_patch(poly)

# set axes
for ax in axes:
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.xaxis.set_ticks_position('bottom')
    ax.set_yticks([])
#  set axes   
axes[0].set_xlabel('x', fontsize=20)
axes[1].set_xlabel('t', fontsize=20)
plt.show()

  从图形上来看，上限的变换，实际上就是坐标轴的变换（时间流逝的速度的变化）。我们常用的坐标是线性坐标，即$x(t)=t$，而实际上，我们也可以取对数坐标$x(t)=log(t)$，甚至任意坐标$x(t)=g(t)$，反映到变上限积分中就是$\int_a^{g(x)} f(t)\mathrm{d}t$，因为变上限积分本质上是一个函数，它的自变量是上限x，而不是积分元t。从函数复合的角度来看，坐标轴的变换就是自变量的复合。
  此外，下面这个结论非常有用：

F′(x)=ddx∫g(x)af(t)dt=f[g(x)]

$$F'(x)=\frac{\mathrm{d}}{\mathrm{d}x}\int_a^{g(x)} f(t)\mathrm{d}t=f[g(x)]$$
4，换元法
  本质上来说，换元法也是坐标变换，因此，在进行换元时，积分元t被替换成u后（$t=\varphi (u),u=\varphi ^{-1}(t)$），积分限也需要进行相应的替换，x变成$\varphi ^{-1}(x)$。
$$\int_{g(a)}^{g(x)} f(t)\mathrm{d}t=\int_{g[\varphi ^{-1}(a)]}^{g[\varphi ^{-1}(x)]} f[\varphi(u)]\varphi'(u)\mathrm{d}u$$
  再看上面的图，设 t=4sin(u),u=arc