用Python学《微积分B》（定积分）

  本文主要学习《微积分B》第7章——“定积分的概念和存在条件”，结合课程中的知识进行一些扩展，并用Python辅助求解课后练习题。关于定积分（Definite Integral）的概念，课本中只介绍了“黎曼积分”（Riemann Integral），而wiki上对定积分的介绍更全面：Integral-Wikipedia。此外，再介绍一个好玩的数学知识讲解网站：Math is fun，它上面对Integration和Definite Integral进行了非常形象的讲解，值得一看。

一、知识点

1，Riemann Integral

  Riemann积分是从求曲线在某区间内与x轴围成的区域的面积出发，分为五个步骤推演出来的。这五个步骤分别是（详见Integral-Wikipedia）：
1）分割
  如上图所示，将区间[a, b]“任意”分割成n份。其中，a=x0<x1<x2<...<xn−1<xn=b。分割是任意的，从上图也可以看出，各分割区域（子区间）并不相等。我们先引入一个描述分割程度的量λ=max{Δxi}，其中，Δxi表示各子区间的大小。
2）取点
  在各子区间上任取一点，如上图所示，有些点取在子区间的左端点（x∗1,x∗5），有些点取在子区间的右端点（x∗2,x∗3），有些取在任意中间点（x∗4）。需要注意的是：Riemann积分的推演，在分割和取点时，都是任意的。网上很多图片都没有真正反映这两个“任意性”，要么分割时进行的是等分，要么取点只取左端点或右端点。这张图是我难得才找到的，能够很好反映这两个“任意性”的图。
3）近似
  在计算各子区间所围成的面积时，将它近似地当作一个{height:f(ξi),width:Δxi}的矩形。这就是“矩形逼近”法。为什么能够这样近似，后面我将给出一个“粗略的证明”。
4）取极限，求和

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

  取极限，求黎曼和（Riemann Summation），这就是Riemann积分的定义。事实上，积分符∫就是来自“Summation”的首字母。
5）证明“极限值与分割的任意性和取点的任意性无关”
  这个证明其实非常重要，只有保证“这个极限值与分割的任意性和取点的任意性无关”，才能说“这个极限值等于函数f(x)在区间[a, b]的定积分”。当然，要直接去证明“任意性”是很难的，因此，后面将给出一个“可积”的定理。


2，“矩形逼近”
  现在，我来粗略地证明“为什么可以用矩形来近似各子区间围成的面积”，即证明“矩形逼近”的合理性。这个证明的过程，我把它分为两步：
1）“梯形近似”
  在前面讲“函数微分的定义”（The definition of ‘Differential’，详见《微积分B》4-2节）的时候，提到：曲线f(x)在x0的邻域内可以用点(x0,f(x0))的切线来近似，（Linear Approximation），事实上，这就是“一阶Taylor公式”近似。
  如果我们把每个子区间的曲线段都用它在该区间上任意一点的切线来近似替代，每个子区间所在区域就变成了“梯形”（Trapeziod），这就是“梯形近似”。更通过更精确地描述来计算这种近似的误差。

Δy=f′(x)Δx+O[(Δx)2]

  其中误差（error）ϵ=O[(Δx)2]。而每个子区间的面积ΔSi与Δy是一次线性关系，即ΔSi∼Δy，因此，梯形近似后，个子区域面积的误差ϵi=ϵ=O[(Δx)2]。
  而总面积S=∑Si，那么总面积的误差

ϵS=n∗ϵi=1Δx∗O[(Δx)2]=O(Δx)

  注：这里为了方便计算，假设子区域误差是相等的。
  所以，当Δx→0时，ϵS→0，即误差趋于0。
2）矩形近似等于梯形近似
  梯形面积计算公式和矩形面积公式分别是：

Sti=a+b2∗Δxi,Sri=f(ξi)∗Δxi

为了简化计算，我只取左端点，计算两种近似之间的误差如下：

Sti=a+b2∗Δxi=[f(xi−1)+dyi2]∗Δxi=[f(xi−1)+f′(ξi)2Δxi]∗Δxi=f(xi−1)∗Δxi+O[(Δxi)2]

显然，各子区域在两种近似之间的误差ϵi=Sti−Sri=O[(Δx)2]，同样的，总面积误差为：

ϵS=n∗ϵi=1Δx∗O[(Δx)2]=O(Δx)

  综上，矩形近似等于梯形近似。


3，“可积”的定理与推论
1）不可积的反例
  “Dirichlet Function”在区间[0, 1]上的积分：

D(x)={10x∈Qx∈R∖Q

  此时，有两种取点方式：全取有理点或全取无理点，如下：

limx→0+D(ξi)∗Δxi{ξi∈[xi−1,xi]=1ξi∈[xi−1,xi]=0∈Q∈R∖Q

2）可积的等价定理
  f为[a, b]上的有界函数，

f∈R[a,b]⇔limλ→0+∑i=1nΔωi(f)Δxi=0

其中Δωi(f)称为振幅。
3）推论
  C[a,b]⊂R[a,b]
  f在[a, b]上只有有限个间断点，则f∈R[a,b]
  [a, b]上的单调函数可积。
4，不定积分与定积分的异同
  不定积分与定积分在形式上（表达式和计算方法）是一致的，但具体的物理意义却不同。
  不定积分一般看作求导的逆运算，“Math is fun”上用“水龙头与水缸”（Tap and Tank）的例子很好地说明了积分与求导之间的关系。它甚至解释了原函数族中的“C”的物理意义（水缸中的存量水）。
  定积分源于计算“曲线在某区间上与x轴围成的面积”。
5，微分与定积分的几何意义
  微分就是在小范围内用切线近似曲线（Linear Approximation），这是一种降维思想。导函数则是各点切线斜率与x轴坐标的映射关系。
  定积分源于计算“曲线在某区间上与x轴围成的面积”，它是面积与x轴坐标区间的映射关系。
  知乎上有人对“微分”与“积分”进行了通俗的比较，可以参考。
6，定积分的性质
  定积分的性质总结为三大类共9个性质，可以参考wiki。所有的这些性质，从面积法的角度来理解，还是很容易的。
1）线性关系：被积函数的加法、数乘的组合。
2）不等式关系：有界性、保号性、保序性等。
3）积分区域的性质：加法、反号等。

二、思考题

1，判断下列两个函数在[-1,1]区间是否可积？

f(x)={x+1,ex,x≤0x>0

g(x)={x,ex,x≤0x>0

答：本题主要根据“可积”的等价定理及其推论来解答
  f(x)在 x = 0 这一点连续，limx→0−f(x)=1=f(0)=limx→0+f(x)。根据推论一，连续函一定可积。
  g(x)除 x = 0 这一点外，都是连续的，且 x = 0 这一点是g(x)的跳跃间断的。根据推论二，它也是可积的。此外，需要注意的是，原题中g(x)的表达式是有误的， x = 0 时，g(x)既等于0也等于1，违反了函数映射关系。
  下面从图形（面积）上来直观的理解g(x)为什么可积：

%matplotlib inlineimport numpy as npimport matplotlib.pyplot as pltfrom matplotlib.patches import Polygona, b, c = -1, 0, 1    # integral limitsx1 = np.linspace(-2, 0)y1 = x1x2 = np.linspace(0, 2)y2 = np.exp(x2)fig, ax = plt.subplots()plt.plot(x1, y1, 'b', linewidth=2, zorder=10)plt.plot(x2, y2, 'g', linewidth=2, zorder=10)#plt.grid(True, zorder=5)# make the shaded regionix1 = np.linspace(a, b)iy1 = ix1verts = [(a, 0)] + list(zip(ix1, iy1)) + [(b, 0)]poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')ax.add_patch(poly)ix2 = np.linspace(b, c)iy2 = np.exp(ix2)verts = [(b, 0)] + list(zip(ix2, iy2)) + [(c, 0)]poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')ax.add_patch(poly)l = plt.axhline(color='red')#plt.text(0.5 * (a + b), 30, r"$\int_a^b f(x)\mathrm{d}x$",#         horizontalalignment='center', fontsize=20)plt.figtext(0.9, 0.05, '$x$')plt.figtext(0.1, 0.9, '$y$')ax.spines['right'].set_visible(False)ax.spines['top'].set_visible(False)ax.xaxis.set_ticks_position('bottom')ax.set_xticks((a, b, c))ax.set_xticklabels(('$-1$', '$0$', '$1$'))ax.set_yticks((0,))ax.set_yticklabels(('$0$',))plt.show()

2，证明：

limn→+∞∫10xnx2+1dx=0

证明：根据微分中值定理可得

∫10xnx2+1dx=ξnξ2+1,ξ∈(0,1)

所以

limn→+∞∫10xnx2+1dx=limn→+∞ξnξ2+1=0

三、选择题

  下面先用python直接给出各个求定积分的题的答案。

from sympy import *init_printing()#Exercise 7-1-1x = Symbol('x')y = integrate(x * atan(x), (x, sqrt(3) / 3, sqrt(3)) )y, y > pi / 9, y < 2 * pi / 3

(−3√3+5π9,True,True)

#Exercise 7-1-2x = Symbol('x')n = Symbol('n', integer=True)f = x ** n / (x ** 2 + 1)limit(integrate(f, (x, 0, 1 / 2)), n, +oo)

0

Exercise 7-1-3

limn→+∞∫10xnexex+1dx=?

解：这道题直接用sympy解不出来，可以参考思考题2，用微分中值定理解答。
Exercise 7-1-4
设f(x)为(0,+∞)上的单调减函数，比较下面几个式子的大小：

∫n1f(x)dx,∫n+11f(x)dx,f(1)+∫n1f(x)dx,∑k=1nf(k)

解：本题主要从面积法来分析。
  先看1式和2式，从面积法很容易理解，它们之间的大小关键看f(x)是否穿过x轴。
  对于最后一个式子，可以将f(k)看作“阶梯函数”，将求和转变为积分。然后，再应用面积法，将对积分的比较转变为对函数曲线的比较。于是，有如下关系：

∫n+11f(x)dx≤∑k=1nf(k)≤f(1)+∫n1f(x)dx

Exercise 7-1-5

limn→+∞∑k=1n2nn2+k2=?

解：这一题直接用sympy的summation函数也是解不出来的。同样地，我们需要利用“阶梯函数”将求和转变为积分。
  函数f(x)=nn2+x2显然是一个单减函数，可以直接利用上一题的结论。

∫n+11nn2+x2dx≤∑k=1nnn2+k2≤f(1)+∫n1nn2+x2dx

  又∫nn2+x2dx=arctan(xn)+C，根据“Newton-Leibniz公式”，可得

arcarctan(xn)∣n+11≤∑k=1nnn2+k2≤f(1)+arcarctan(xn)∣n1

  然后再用极限的“夹逼定理”，两边极限值都等于arctan(1)=π4
解2：利用极限的定义，进行不等式缩放

12=∑1nnn2+n2≤∑k=1nnn2+k2≤∑1nnn2+12≤∑1nnn2=1

  则S(n)=∑nk=1nn2+k2是单调有界数列，必定收敛。（只能到这一步了，待高手）
  下面结合这两题，用图形法来直观理解：

import numpy as npimport matplotlib.pyplot as pltimport matplotlib.ticker as tickerdef func(x):    n = 10    return n / (n ** 2 + x ** 2)a, b = 2, 9  # integral limitsx = np.linspace(0, 11)y = func(x)x2 = np.linspace(1, 12)y2 = func(x2-1)ix = np.linspace(1, 10, 10)iy = func(ix)fig, ax = plt.subplots()plt.plot(x, y, 'r', linewidth=2, zorder=15)plt.plot(x2, y2, 'g', linewidth=2, zorder=15)plt.bar(ix, iy, width=0.95, color='gray', align='edge', ec='olive', ls='--', lw=2,zorder=10)plt.ylim(ymin=0)plt.figtext(0.9, 0.05, '$x$')plt.figtext(0.1, 0.9, '$y$')ax.spines['right'].set_visible(False)ax.spines['top'].set_visible(False)ax.xaxis.set_major_locator(ticker.IndexLocator(base=1, offset=1))plt.show()

  如上图所示：红色曲线表示函数f(x)，它在(0,+∞)是单减函数，它在[1, n]区间上的积分∫n1f(x)dx就是红色曲线在区间[1, n]之间与x轴围成的面积；∑nk=1f(k)表示对应的矩形面积和。很明显，矩形面积和大于或等于曲线在对应区间与x轴围成的面积。也许你会说，明明只看到大于，怎么会等于呢？那是因为我画图时，取的 n = 10，随着n取值的增大，超出的面积会越来越小，当n→+∞时，是有可能相等的。
  再看上图的绿色曲线，它表示函数f(x-1)。很明显，矩形面积和是小于或等于绿色曲线的积分的。下面通过公式推演，将绿色曲线的积分与最右边的式子联系起来：

∫n+11f(x−1)dx=∫n0f(t)dt=∫n1f(t)dt+∫10f(t)dt

  现在，只需证明下式即可

limn→+∞∫10f(t)dtf(1)=1

  应用微分中值定理：

limn→+∞∫10f(t)dtf(1)=limn→+∞f(ξ)f(1)=C,ξ∈(0,1)

  对于f(x)=nn2+x2，是很容易证明的这个极限值为1。
  另解：前面我们在讲“Riemann Integral”时，粗略地证明了“用矩形来逼近积分面积”，其中矩形的取法可以是左端点，也可以是右端点，即“取点的任意性”。事实上，这就是“Darboux Integral”。不论怎样的取法，不影响极限值（积分结果），从图形上来看，相当于保持曲线不动，左右移动矩形（左端点、右端点），极限值相等。而这道题刚好是反过来，相当于保持矩形不动，左右移动曲线，同理可得它的极限值也是相等的。


Exercise 7-1-6

∫1−1[sin5(x)+1+x2−−−−−√]dx=?

解：根据定积分的线性组合性质，可以把这个定积分分为“三角函数部分积分”+“无理式部分积分”，其中，三角函数部分是奇函数，在对称区间的积分，等于0，而无理式部分的积分是一个简单二次根号，很容易想到“三角函数平方公式”进行有理化，如下：

∫1−11+x2−−−−−√dx=∫π2−π2cos(θ)d[sin(θ)]=∫π2−π2cos2(θ)dθ=π2

Exercise 7-1-7

∫1−1x2ln(x+1+x2−−−−−√)dx=?

解：先分部积分，再三角函数平方公式有理化

∫1−1x2ln(x+1+x2−−−−−√)dx=13[x3∗ln(x+1+x2−−−−−√)∣1−1−∫1−1x3d[ln(x+1+x2−−−−−√)]=13[x3∗ln(x+1+x2−−−−−√)∣1−1−∫1−1x3x2+1−−−−−√dx=13[x3∗ln(x+1+x2−−−−−√)∣1−1−∫1−1x3x2+1−−−−−√dx]

  前一部分计算化简后等于0，后一部分是“奇函数在对称区间积分等于0”


Exercise 7-1-8
比较两个积分的大小，直接画图如下：

import numpy as npimport matplotlib.pyplot as pltx = np.linspace(0, np.pi / 2)y = np.sin(np.sin(x))y2 = np.cos(np.sin(x))plt.plot(x, y, 'r', zorder=10)plt.plot(x, y2, 'g', zorder=10)plt.grid(True, zorder=5)plt.show()