用Python学《微积分B》（微积分应用）
微积分是现代科学中一种核心工具和基本概念，在多个领域发挥着重要作用。例如，在几何学中用于计算平面面积与曲线长度，在立体几何中则用于确定空间图形的体积与旋转曲面面积。此外，在物理学研究中也广泛应用"微元法"进行建模与计算。因此掌握微积分的核心要诀，则是本章学习的关键所在。深入理解其背后的理论体系至关重要：首先需要透彻掌握"导数-Derivative"与"不定积Integral"的基本定义与内涵；其次要认识到这一理论体系在英语中有专门的称谓："Infinitesimal calculus"亦即"无穷小量演算"。从这个角度而言它不仅是一种数学工具更是一个揭示微观规律并最终实现宏观计算的思想体系。具体来说可以通过以下步骤进行学习：首先将复杂问题进行分割其次选取恰当的研究对象对局部现象进行深入分析最后将微观规律推广至宏观场景并评估误差范围以确保结果精确度。
关于定积分的应用部分建议参考以下资源获取详细讲解：
AREAS VOLUMES OF SOLIDS OF REVOLUTION LENGTH OF A CURVE AREAS OF SURFACES OF REVOLUTION WORK FLUID PRESSURE
以及
Integral Application
这些材料都提供了非常详尽且易于理解的内容值得深入研究。
一、分割
分割是微积分方法的第一步，也是微积分应用中非常重要的一步。算法中有“分而治之”的策略（Divide-and-conquer algorithms），微积分的“分割”也正暗合这种思想。另外所谓“微观化”，通俗理解就是取待研究的对象的一小部分作为单元，放大了仔细研究，找出特征，然后再总结整体规律。而微积分的“分割”也正是这个“取一小部分作为单元”。
普遍来说，有两种分割方式：直角坐标系分割和极坐标系分割。
1，直角坐标系分割
对于直角坐标系分割，我们已经和熟悉了，前面将定积分定义的时候，就是在直角坐标系下用“矩形逼近”的方法来计算曲线与x轴围成的面积。它是沿x轴分割成n小段\{\Delta x_i \}，即在直角坐标系下分割是按自变量进行分割。
当然，直角坐标系下也可以沿y轴分割，本质上，直角坐标系中沿x轴分割和沿y轴分割意义是一样的。将沿y轴分割看作是：x=f^{-1}(y)，将函数关系反转，同时也将坐标轴反转。
2，极坐标系分割
同样地，极坐标也是按自变量分割。只是，直观上看，与直角坐标系的分割差异较大。如下图：
总结：
无论采用何种坐标系，在进行单元划分时都是按自变量进行操作。这种划分是由函数映射关系所决定的；一旦选定合适的坐标系后，在已知自变量的情况下可以通过函数运算得到对应的函数值。从图形上看，在这种划分下每个分割单元中‘不规则边’的数量最少，并且最好是只有一条不规则边。因此，在研究实际问题建模的时候，“选择合适的一种坐标系”是关键
二、取点
由积分的定义可知，在选取区间内任意一点时，并不影响积分的结果。然而，在实际应用过程中为了便于计算或进行定性分析等目的，在选取这些特定点时通常会选择一些典型位置的点值来代表整个区间上的函数值变化情况。例如，在证明某个不等式成立的过程中, 我们会将左边和右边端部对应的表达式分别视为两条曲线下的积分, 而中间求和项对应的矩形面积之和则与这两条曲线围成区域相对应, 这样一来每个矩形的一对相邻端点就会分别位于这两条不同的曲线上。
此外，Darboux Integral也是取得左右两端点。
三、近似
近似是微积分方法最重要的一步。通过“分割”，有了微观上的“单元”后，这个“单元”还是不太适合直接研究，因为它不规则，只有通过近似，将这个不规则的“单元”近似为一个“规则的单元”，这样才能继续下一步研究。这么说来，“近似”是整个微积分中最有创意，最需要发挥人的联想能力的一步。
1，不规则近似为规则
可以这么说：近似就是在微观上将不规则的“单元”替换为规则的“单元” 。回到面积法，我们无法直接计算一个曲边图形的面积，但是在微观单元上，我们可以用一个相似的直边图形来替代它。直观上看，只要这个微观单元足够小，这个替代的误差也就足够小。也就是说，这个替代某种意义上是可行的，误差是可控制的。前面在讲“分割”的时候说“要使不规则的边的数量尽可能小”，实际上也就是方便做“近似”。
2，直线代曲线
更具体一点，在前面介绍的定积分定义和上面的“极坐标图形”问题，可以发现这两类问题近似的实际就是“用直线替代曲线”。在直角坐标系中，这么替代后，分割单元由曲边梯形变成了斜边梯形；极坐标系中，这么替代后分割单元由曲边三角形变成了普通三角形。这一步做完后，你会发现在微观上，原来不可计算的问题变成了可计算的问题了。
注意，在极坐标系中，计算面积时，既可以用“三角形近似”（Triangle），也可以用“圆弧近似”（Arc），后面将讨论这两种近似的误差是一致的。
3，套用计算公式
之所以要将不规则单元替换为规则单元，是因为规则单元可以套用计算公式。
替换完成后，下一步就是针对待求解的问题，对“规则单元”套用已知的公式。待求解的问题不同，套用的公式显然也不同。比如：
1）Riemann和定义的例子
待求解的是在区间[a,b]上曲线与x轴围成的面积，因此套用的是平面面积公式。
2）该区域的极坐标曲线积分（上图）。
待求解的是在区间[\theta_1, \theta_2]上曲线所围区域的面积，并应用圆弧区域面积计算公式进行计算。
3）平面曲线弧长
平面曲线在微观层次上近似表现为一系列微小的折线段（即切线），其计算方法基于勾股定理的原理：将曲线分割成无数个极小的直角三角形，并通过计算这些三角形斜边长度之和来确定整体弧长。也就是说，在极限情况下这一过程趋近于精确值。
4）极坐标曲线长度
请务必注意的是，在推导过程中涉及到了\pi这一常数，并且该公式本身是通过近似方法得到的。同样地，请记住我们也可以将其推广至旋转体的体积与表面积计算。
四、求和
前几步均处于微观层面上，在这种情况下仅当采用Riemann和时才能返回到宏观层面
其中，F_i表示各个微观单元的公式
五、误差分析
“近似”是发挥人联想能力的时候，但联想完了之后，我们要证明这种“近似”是可行的，即证明“误差在可接受范围内”。当然，对于误差的计算是要回到宏观层面上来的。一是我们原本要研究的就是一个宏观问题，最后的计算结果只有回到宏观上看才有意义；二是微观上的小误差有可能累积到宏观上变成大误差，正所谓“差之毫厘谬以千里”。
1，平面曲线积分误差分析
在“定积分”那一节，我通过“无穷小的运算”证明了“梯形近似”的误差 \epsilon=O(\Delta x) ，同时也证明了“矩形近似与梯形近似的误差在同一个级别—— O(\Delta x) ”。传送门-Integral
2，极坐标曲线积分误差分析
现在我们来证明极坐标曲线积分的“三角形近似”和“圆弧近似”的误差 。
1）圆弧近似（Arc）
先用弧线代曲线
误差为 \epsilon_l = O[(\Delta \theta_i)^2]
再算单元面积
根据单元面积公式可知，
S_{i}与\Delta\theta_{i}呈一次线性关系，
即有 S_{i}\sim\Delta\theta_{i}。
那么用弧形面积近似后产生的误差为 \epsilon_{i}=O[(\Delta\theta_{i})^{2}]
求和后的总误差为
在计算过程中，默认各子区域的误差相等。
由此可知，在Δθ趋近于零的情况下，
\epsilon \rightarrow 0
第2节 三角形近似（Triangle）
采用直线逼近曲线的方式重新计算面积
同样地，可以计算误差为：
根据无穷小的替换
此外，通过无穷小的替换，也可以证明这两个面积相等（Riemann和）
关于三角形和圆弧逼近的关系论述可能存在一定的因果关系颠倒情况；不过这种讨论有助于我们更好地理解相关概念。实际上，在本课程的第一章中介绍割圆术时……扈老师通过极限的方法演示了两者在极限条件下相等性原理；如果有兴趣可以回看相关教学视频以进一步加深理解。
如下图，直观上都可以感受到这个误差之大。
当我们计算曲线的长度时，在某些情况下我们需要采用特定符号替代整个曲线。具体而言，并基于毕达哥拉斯定理（Pythagoras theorem），我们将该符号转换为... 如下所示：
2）“”代曲线的误差计算
每个单元误差为 \epsilon_i = O[(\Delta l_i)^2]
求和后总误差为
为了解决问题，在分析中我们假设各小段的误差相等。
随着Δl趋近于零时，
3）两种近似的比较
六、课后练习
Exercise 7-5-1
求由双纽线 \rho^2 = 2a^2\cos(2\theta) 所界定的平面区域的面积
解：观察到该曲线具有对称性特征，在其整体结构中选取一个对称的部分进行详细推导即可确定整个区域的面积值
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
alpha = 1
theta = np.linspace(0, 2*np.pi, num=1000)
x = alpha * np.sqrt(2) * np.cos(theta) / (np.sin(theta)**2 + 1)
y = alpha * np.sqrt(2) * np.cos(theta) * np.sin(theta) / (np.sin(theta)**2 + 1)
plt.plot(x, y)
plt.grid()
plt.show()
#Exercise 7-5-1
from sympy import *
init_printing()
t, a = symbols('t a')
f = a ** 2 * cos(2 * t)
4 * integrate(f, (t, 0, pi / 4))Exercise 7-5-2
计算由心形线（heart-line） \rho = a * (1 +cos(\theta)) 所围成的平面区域的面积
解答：观察到心形线是一个对称图形，在其基础上只需计算其一半区域即可得出总面积值。需要注意的是，在实际绘图过程中虽然图形可能不完全与公式对应，但按照给定公式绘制的心形曲线通常会呈现出水平偏右的特点。
import matplotlib.pyplot as plt
import numpy as np
t = np.arange(0,2*np.pi, 0.1)
x = 16*np.sin(t)**3
y = 13*np.cos(t)-5*np.cos(2*t)-2*np.cos(3*t)-np.cos(4*t)
plt.plot(x,y)
plt.grid()
plt.show()
#Exercise 7-5-2
from sympy import *
init_printing()
t, a = symbols('t a')
rho = a * (1 + cos(t))
integrate(1 / 2 * rho ** 2, (t, 0, pi)) #Exercise 7-5-3
from sympy import *
init_printing()
x = Symbol('x')
f = x *
g = sqrt(x)
a, b = solve(Eq(f, g), x)
a, b, integrate(f - g, (x, a, b)) #Exercise 7-5-4
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
a = 1
t = np.linspace(0, 2 * np.pi)
x = a * (t - np.sin(t))
y = a * (1 - np.cos(t))
plt.plot(x, y)
plt.grid()
plt.show()
Exercise 7-5-4
求旋轮线的弧长
解答：虽然题目看似简单, 但实际操作中发现无法通过sympy直接求解该曲线的弧长。
根据弧长积分公式\mathrm{d}s = \sqrt{1 + \left( \frac{\mathrm{d} y}{\mathrm{d} x} \right)^2 } \,\mathrm{d}x
将积分元替换为t得：
再次利用三角函数关系，进行有理化
#Exercise 7-5-5
from sympy import *
init_printing()
t, a = symbols('t a')
r = a * (1 + cos(t))
f = sqrt(r ** 2 + diff(r, t) ** 2)
f需要注意的是，在处理函数f时，无法直接通过Sympy进行积分运算；类似于上一题的做法进行化简处理后可得积分为8a
#Exercise 7-5-6
from sympy import *
init_printing()
t, a, c = symbols('t a c')
x = a * cos(t)
y = a * sin(t)
z = c * t
f = sqrt(x.diff(t) ** 2 + y.diff(t) ** 2 + z.diff(t) ** 2)
integrate(f, (t, 0, 2 * pi)) #Exercise 7-5-9
x = Symbol('x')
f = pi * sin(x) *
integrate(f, (x, 0, pi)) #Exercise 7-5-10
x = Symbol('x')
y = x ** 3 / 3
perimeter = 2 * pi * y
Dlx = sqrt(1 + diff(y, x) ** 2)
f = perimeter * Dlx
integrate(f, (x, 0, 1))Exercise 7-5-11
计算由平面曲线图形 \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} = 1所围区域的面积。
解答过程如下：该曲线属于椭圆图形。其中a、b分别代表长半轴与短半轴长度。采用参数方程法更为简便。
#Exercise 7-5-11
from sympy import *
init_printing()
t, a, b = symbols('t a b')
x = a * cos(t)
y = b * sin(t)
f = y * diff(x, t)
integrate(f, (t, 0, 2 * pi)) #Exercise 7-5-12
from sympy import *
init_printing()
x = Symbol('x')
y = 1 - x *
a, b = solve(y, x)
f = pi * y *
integrate(f, (x, a, b)) #Exercise 7-5-13
from sympy import *
init_printing()
x = Symbol('x')
y = sin(x)
f = pi * y *
integrate(f, (x, 0, pi)) #Exercise 7-5-14
from sympy import *
init_printing()
y = Symbol('y')
eq = Eq(y ** 2 / 2, y + 4)
a, b = solve(eq)
f = y ** 2 / 2 - (y + 4)
a, b, integrate(f, (y, a, b)) %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
y = np.linspace(-3, 5)
x1 = y ** 2 / 2
x2 = y + 4
plt.plot(x1, y)
plt.plot(x2, y)
plt.grid()
plt.show()