【分　 类】 基础科学
【关 键 词】
【来 　源】 互联网
【收　 录】

正文：

冯长建，孙禹辉，吴鑫辉

（大连民族大学机电工程学院，辽宁省大连市 116600）

摘要：第一类和第二类曲面积分分别称为标量函数的曲面积分和矢量场的曲面积分，两者既有区别，也有联系。这里从不同形式表达的曲面的面积微元入手，推导出适用于两类曲面积分的计算公式。针对不同类型的曲面，利用Python Sympy（符号计算库）进行了数学建模实验，实验结果表明这些算法易于计算机代码实现，且计算准确、高效。

关键词： 矢量场；曲面积分；重积分；数学实验；Python

中图分类号：O172 　　　　　　　文献标志码：A

　　

　　两类曲面积分概念抽象，计算复杂，是高等数学的教学难点所在。随着计算技术的迅猛发展，新工科时代的高等数学也要与时俱进，要能够适应计算机运算、符号理解、甚至公式推导，从而降低人工计算的繁杂性。蓬勃发展计算机符号系统，也称为计算机代数系统，已经可以解决高等数学中的绝大部分问题。这里从曲面的矢量方程出发，导出了适宜于计算机计算的两类曲面积分的公式，并利用Python中Sympy（符号计算库）进行数学实验验证。实验表明，算法精确、有效。Python Sympy模块虽然还不够成熟，但它是免费、开源的系统，研究它在高等数学中的应用，对进一步推动国内计算机代数系统的开发具有积极意义。

1 两类曲面积分计算原理

　　如果是标量函数，则第一类曲面积分，又称为对曲面面积的曲面积分，或标量函数的曲面积分，计算公式定义[1,2,3]如下，

　　. （1）

　　如果是定义在三维空间上的矢量场，它的各个分量定义在光滑曲面上并连续，选择曲面的单位法向量来定向曲面。那么在整个曲面的曲面积分定义[1,2,3]为，

（2）

该积分又称为矢量场穿越曲面的通量，称为第二类曲面积分，矢量场的曲面积分，或对坐标的曲面积分。表示曲面的面积微元。也可以认为是直接在矢量曲面上积分。

1.1 参数化曲面的曲面积分公式

　　如果一个曲面由参数化区域上的矢量方程[1,2,3]给出，

　　 （3）

则它的面积微元[1]由下式给出，

　　, （4）

　　曲面定向的单位法向量的表达式由下式计算，

　　, （5）

　　根据第一类曲面积分的定义公式（1），则第一类曲面积分可如下计算，

. （6）

　　根据第二类曲面积分的定义公式（2），第二类曲面积分计算公式如下，

. （7）

公式中的为曲面的投影平面区域。

　　如果发现曲面的定向和规定的不一致，只要交换向量积的顺序即可。一般而言，曲面法向量指向外侧则曲面积分取正号，反之取负号。

1.2 显式曲面的曲面积分公式

　　对于由给出的显式曲面，利用参数化曲面的定义，对于曲面定义在区域上，而函数具有连续偏导数，如果选择和作为参数，则参数方程可以写为：

　　　 （8）

那么，

（9）

（10）

　　曲面的面积微元表达为，

（11）

　　于是显式曲面的第一类曲面积分的计算公式（1）可以写为，

　　 （12）

　　于是第二类曲面积分的计算公式（2）可以写为，

　　 （13）

这里的为曲面在面上的投影区域，而为其面积微元。

　　实际计算上，可以构造一个三元函数，则其梯度为，

　　 （14）

　　于是，公式（12）可以简单写为，

　　 （15）

　　注意，第二类曲面积分，曲面是有方向的，曲面的法线方向和坐标轴方向一致为曲面积分为正，否则取符号。也就是曲面的上侧、前侧、右侧取正。

1.3 隐式曲面的曲面积分公式

　　　比如对于常数，方程所表达的曲面称之为隐函数定义的曲面，不妨简称为隐式曲面。定义矢量为垂直于隐式曲面下方的平面区域的单位法向量。假定曲面光滑（函数可微，曲面的梯度非零，且连续），且。

　　假定矢量为单位向量，于是方程所确定的曲面的平面区域就在面上。根据假定，在曲面上有。隐函数理论暗示着曲面的图形就是可微函数的图形，不过这里的并不能显式已知。定义参数和，并令，那么。于是，

　　　 （16）

给出了曲面的参数化方程。可以根据参数化曲面的面积微元计算方法进行计算。首先计算矢量函数的偏导数，得到，

（17）

　　利用隐函数的求导法则，可以获得求得偏导数，

　　 （18）

　　把这些偏导数代入到公式（16）中可得

　　　 （19）

　　利用矢量积运算可得，

　　　 （20）

　　于是，隐函数曲面的面积微元可由下式给出，

. （21）

　　法线单位法向量可表示为曲面函数G梯度的归一化向量，

　　. （22）

　　于是隐式函数确定的曲面的第一类曲面积分的计算公式为，

　　 （23）

　　第二类曲面积分的计算公式为，

　　 （24）

　　这里第二类曲面积分的正符号规定和1.2节的一致。

2 二重积分和区域积分的Python实现

　　　两类曲面积分传统计算方法，都是转化为平面区域上的二重积分进行计算。Python Sympy计算二重积分主要有两种方法，第一种是把二重积分转化为二次累次积分，第二种方法上采用矢量积分，直接在区域上积分[4]。除此之外，Python可以在曲面上直接建立曲面积分，这样一来，可以大大简化计算步骤，而不必转化为平面区域上的二重积分。实践上发现直接曲面积分，算法更加稳定，其核心任务就是能够恰当地描述曲面的参数化区域。下面举例说明Python中的符号积分计算方法。

　　　例1 计算，其中

　　　解：下面利用多种方法求解此二重积分问题。

　　　解法1：直接在直角坐标系下累次积分，Python实现代码如下。

　　　In[1]: from sympy import *#调用符号计算库

　　　x,y=symbols('x,y') #声明符号变量

　　　R=symbols('R',positive=True) #声明符号常量为正

　　　intfun=x**2+y**2#被积函数表达式

　　　integrate(intfun,(y,-sqrt(R**2-x**2),sqrt(R**2-x**2)),(x,-R,R))

　　　#累次积分

　　　Out[1]

　　　解法2：利用极坐标转化积分区域，再进行二次积分。

　　　In[2]: from sympy import *#调入符号运算库

　　　x,y,t=symbols('x,y,t') #定义符号变量

　　　r,R=symbols('r,R',positive=True) #定义正的符号变量

　　　intfun=x**2+y**2#被积函数表达式

　　　intfun=intfun.subs({x:r*cos(t),y:r*sin(t)}) #极坐标替换

　　　intfun=intfun.simplify()#尝试化简

　　　integrate(intfun*r,(r,0,R),(t,0,2*pi)) #二次积分

　　　Out[2]

　　　解法3： 利用参数化描述平面积分区域，直接进行区域积分

　　　In[3]: from sympy import *#调入符号运算库

　　　from sympy.vector import *#调入矢量模块

　　　C=CoordSys3D('C') #定义坐标系

　　　x,y,t=symbols('x,y,t') #定义符号变量

　　　r,R=symbols('r,R',positive=True) #定义正的符号变量

　　　x=r*cos(t) ;y=r*sin(t) #曲面参数方程的x和y分量

　　　intfun=C.x**2+C.y**2#被积函数

　　　Dxy=ParametricRegion((x,y),(r,0,R),(t,0,2*pi)) #圆域的参数化表达

　　　vector_integrate(intfun,Dxy) #矢量积分

　　　Out[3]

　　　In[]和Out[]分别表示程序的输入和输出，其中中括号中的数字表示序号，以下类似。从计算效果看几种方法输出结果一致，且都是精确的符号解。

　　　例2 第一类曲面积分计算，其中为上半球面。

　　　解：通过该例子，说明曲面区域直接积分的使用方法。

　　　In[4]: C=CoordSys3D('C') #定义坐标系

　　　x,y,z,u,v=symbols('x,y,z,u,v') #定义符号变量

　　　H=C.z**3#被积函数

　　　x=sin(u)*cos(v)#定义球坐标系x分量

　　　y=sin(u)*sin(v) #定义球坐标系y分量

　　　z=cos(u) #定义球坐标系z分量

　　　semisphere=ParametricRegion((x,y,z),(u,0,pi/2),(v,0,2*pi)) #半球区域参数化表达

　　　vector_integrate(H,semisphere) #半球面执行标量曲面积分

　　　Out[4]

　　　例3 第二类曲面积分计算，其中Σ为曲面在第1卦限的部分取上侧。

　　　解：通过该例子，说明第二类曲面积分也可以直接建立参数化曲面上，Python程序如下。

　　　In[5]: C=CoordSys3D('C') #定义坐标系

　　　x,y,z,t=symbols('x,y,z,t') #定义符号变量

　　　r=symbols('r',positive=True) #定义正符号变量

　　　F=C.x*C.i+C.y*C.j+2*C.z*C.k#定义矢量场

　　　x=r*cos(t) #曲面参数x分量

　　　y=r*sin(t) #曲面参数y分量

　　　z=1-x**2-y**2#曲面参数z分量

　　　surface=ParametricRegion((x,y,z),(r,0,1),(t,0,pi/2)) #抛弃面参数化区域表达

　　　vector_integrate(F,surface) #直接矢量场积分

　　　 Out[5]

　　　例4 计算第二类曲面积分或通量，这里，积分曲面为如图1所示的缺失单位立方体的几何体的整个边界面，方向向外。

图1 例5中的积分区域

　　　解：这里通量如果利用面积微元法不易求解，但由于是封闭曲面，利用高斯定理，可以在两个体积域上积分，然后，求两者之差即可。这里直接建立体积域上的积分。

　　　In[6]: C=CoordSys3D('C')#定义坐标系

　　　x,y,z=symbols('x,y,z') #定义符号变量

　　　F=C.x*C.i+C.y*C.j+C.z*C.k#定义矢量场

　　　cuboid1=ParametricRegion((x,y,z),(x,0,2),(y,0,2),(z,0,2)) #大立方体定义

　　　cuboid2=ParametricRegion((x,y,z),(x,1,2),(y,1,2),(z,1,2)) #小立方体定义

　　　I1=vector_integrate(divergence(F),cuboid1) #散度理论

　　　I2=vector_integrate(divergence(F),cuboid2) #散度理论

　　　Flux=I1-I2#通量计算

　　　Flux#输出结果

　　　Out[6]

3 曲面积分的Python综合实验

　　　通过下面的例题，把上面所述的各种计算公式进行Python编程实验。

　　　例5 计算曲面积分，曲面为抛物面，在面上方的部分，曲面取上侧。

　　　解法1：利用矢量场直接在参数曲面上积分，即利用第1节的公式（2）左侧计算，Python实现代码如下。

　　　In[7]: from sympy import *#调入符号库

　　　from sympy.vector import *#调入矢量库

　　　x,y,z,v=symbols('x,y,z,v') #定义符号变量

　　　u=symbols('u',positive=True) #定义符号变量，并假设为正

　　　F=C.x*C.i+C.y*C.j+C.z*C.k#定义矢量场

　　　G=C.z+C.x**2+C.y**2-1#定义曲面的构造函数

　　　x=u*cos(v)；y=u*sin(v)；z=1-x**2-y**2#定义参数化曲面各分量

　　　surface=ParametricRegion((x,y,z),(u,0,1),(v,0,2*pi)) #参数化曲面

　　　vector_integrate(F,surface) #公式（2）的应用

　　　Out[7]

　　　解法2：把矢量曲面积分转化为第一类曲面积分，利用第1节中的公式（2）的右侧计算，Python实现代码如下。

　　　In[7]: F=Matrix([x,y,z])#定义矢量场向量

　　　gradG=gradient(G) #计算曲面的梯度

　　　n=gradG.normalize()#曲面的单位法向量计算

　　　vector_integrate(F.dot(n),surface) #公式（2）的应用

　　　Out[7]

　　　解法3：利用参数化曲面，即第1.1小节中的公式（7）计算,Python实现代码如下。

　　　In[7]: r=Matrix([x,y,z]) #定义矢量曲面

　　　ru=r.diff(u)；rv=r.diff(v) #计算矢量的偏导数

　　　intfun=F.dot(ru.cross(rv)) #计算被积函数

　　　intfun=intfun.simplify()#化简

　　　integrate(intfun,(u,0,1),(v,0,2*pi)) #二重积分

　　　Out[7]

　　　解法4： 利用显式曲面公式，即1.2小节的公式（15）计算，Python实现代码如下。

　　　In[7]: F=C.x*C.i+C.y*C.j+C.z*C.k#定义矢量场

　　　gradG=gradient(G) #定义构造曲面的梯度

　　　Dxy=ParametricRegion((u*cos(v),u*sin(v)),(u,0,1),(v,0,2*pi)) #平面区域参数化

　　　intfun=F.dot(gradG) #被积函数

　　　intfun=intfun.subs(C.z,1-C.x**2-C.y**2) #变量替换

　　　vector_integrate(intfun,Dxy) #区域积分

　　　Out[7]

　　　解法5：利用隐式曲面公式，第1.3小节中的公式（24）计算，Python实现代码如下。

　　　In[7]: p=C.k#定义定向向量k

　　　intfun=F.dot(gradG/Abs(gradG.dot(p))) #被积函数表达

　　　intfun=intfun.subs(C.z,1-C.x**2-C.y**2) #变量替换

　　　vector_integrate(intfun,Dxy) #积分

　　　Out[7]

　　　当然，实现的方法还有很多，篇幅有限，不再举例，从上述计算结果可以看出，计算结果准确无误。

4 结论

　　从上述的基于Python的数学计算实验结果可以看出，利用CAS（计算机代数系统），编程简单，计算准确、高效。目前，国内的CAS的研究和开发起步较晚，因此研究各种算法的计算机实现途径，尤其是开源系统Python的代码实现，对于解决复杂计算问题具有十分积极的意义，也必将推动我国自主符号计算系统的研发进程。

参考文献:

[1] 宁荣健, 周玲. 曲面积分参数方程计算方法的若干注记[J]. 大学数学,2020,36(5):67-72.

[2] WEIR, HASS, GIORDANO. 托马斯微积分(第11版) (下册) (影印版)[M]. 北京: 高等教育出版社,2016.

[3] 史迪沃特编著. 微积分(第7版)(下册)(影印版)[M]. 北京: 高等教育出版社, 2014.

[4] 周群, 刘雄伟. Mathematica区域构建与多元函数积分计算的直接方法[J]. 教育现代化,2019,6(103):165-168.