多元函数相关问题的Matlab实现

Matlab Implementation of Multivariate Function Related Problems

摘要: 本文借助于Matlab强大的绘图功能以及高效计算功能，针对多元函数偏导数的几何意义以及二重积分、三重积分的计算问题给出了Matlab程序，增强了利用数学软件工具分析问题、解决问题的能力。

Abstract: This paper utilizes Matlab’s powerful drawing and efficient calculation functions to provide Matlab programs for the geometric meaning of partial derivatives of multivariate functions and the calculation problems of double and triple integrals, enhancing the ability to analyze and solve problems using mathematical software tools.

1. 引言

多元函数偏导数以及重积分等知识点是各高校数学类课程的重要内容，是学习其他专业课程的基础。多元函数偏导数在经济学、工程学、流体力学、电磁学以及生物学等多个科学领域都有着重要的应用，比如具有偏导数思想的拉格朗日方法是求解多元函数的最值的一个重要工具，梯度、散度和旋度等都是通过偏导数定义的，同时，偏导数又是构成偏微分方程的重要组成部分；二重积分在几何上用于计算平面区域的面积以及空间体的体积，在物理学中用来计算平面区域的质心与转动惯量、磁通量与能量等；三重积分用于计算空间区域的体积，在物理学中用来确定空间中质量的分布，在电磁学和流体力学中，三重积分用于描述和计算物理场等。但由于多元函数偏导与重积分理论的抽象性，使得学生在学习过程中失去学习兴趣。而Matlab具有强大的科学计算功能，能够提供高精度的数值计算方法，可以准确地进行偏导数、积分等计算；同时又具有强大的绘图功能，可以直观地展示函数图像、积分区域和结果，有助于理解问题的本质，对于理论研究和教学非常有利。

《浙江省教育事业发展“十四五”规划》中指出坚持数字赋能、创新驱动，以促进信息技术与教育教学融合创新为主线，高水平推进教育管理信息化，注重技术应用理论，全面提升师生信息素养[1]。鉴于Matlab软件具有强大的绘图功能以及高效计算功能，可作为数学类课程教学的重要辅助工具[2]-[5]。本文借助于Matlab软件对多元函数偏导数以及重积分知识点进行编程实现，按照先直观后抽象的原则，通过数据的可视化阐明知识点的本质来培养学生应用数学软件解决实际问题的能力。

2. Matlab可视化多元函数偏导数的几何意义

2.1. 预备知识

二元函数$z=f\left(x,y\right)$ 在点$\left(x_0,y_0\right)$ 的偏导数的几何意义如下([6]: p. 66)：

设$M_0\left(x_0,y_0,f\left(x_0,y_0\right)\right)$ 为曲面$z=f\left(x,y\right)$ 上的一点，过$M_0$ 作平面$y=y_0$ ，截此曲面得一曲线，此曲线在平面$y=y_0$ 上的方程为$z=f\left(x,y_0\right)$ ，则偏导数$f_x\left(x_0,y_0\right)$ 就是曲线在点$M_0$ 处的切线对x轴的斜率。同样，偏导数$f_y\left(x_0,y_0\right)$ 的几何意义就是曲面被平面$x=x_0$ 所截得的曲线在点$M_0$ 处的切线对y轴的斜率。

2.2. 二元函数偏导数Matlab程序实现

通过例题，利用Matlab作出函数图像帮助学生理解多元函数偏导数的几何意义。

例1 讨论二元函数$z=2-2x^2-y^2$ 在点$\left(1/2,1/2\right)$ 处的对自变量x的偏导数的几何意义。

分析 首先，我们知道二元函数$z=2-2x^2-y^2$ 表示椭圆抛物面；

其次，在点$\left(1/2,1/2\right)$ 处作平面$y=1/2$ ，与椭圆抛物面$z=2-2x^2-y^2$ 相交形成的曲线为$z=2-2x^2-{\left(1/2\right)}^2=1.75-2x^2$ ；

最后，在$x=1/2$ 处作该曲线的切线，经过计算得到切线方程为$z=2.25-2x$ 。

代码如下：

z1 = @(x,y)x.^2*(-2)-y.^2+2;

x = -2:.05:2; y = -2:.05:2;

[X,Y] = meshgrid(x,y);

Z1 = z1(X,Y);

surf(X,Y,Z1); %绘制椭圆抛物面

hold on;

z = -5:.1:5;

[X,Z2] = meshgrid(x,z);

Y= 0.5+0*X+0*Z2;

surf(X,Y,Z2); %绘制平面

hold on;

X = -2:.05:2;

Z = 1.75 - 2*X.^2;

plot3(X,Y,Z,'LineWidth',2,'Color',[0,0,0]); %绘制相交的曲线

scatter3(0.5,0.5,1.25,40,'filled'); %绘制切点

hold on;

Z4 = 2.25 - 2.*X;

plot3(X,Y,Z4,'LineWidth',1,'Color',[0,1,0]); %绘制切线

text(-1.24,3.15,-13.77,'Z = 1.75 - 2*X.^2 \rightarrow ')

text(1.02,-0.76,8.6,'\leftarrow y =0.5')

text(1.5,-0.8,-3.0,'\leftarrow Z =2-2*X.^2-y^2')

text(1.72,0.03,-0.22,'\leftarrow z=2.25-2*x')

输出函数图像如下：

Figure 1. Image of geometric meaning of partial derivative in Example 1

图1. 例1中偏导数几何意义图像

从图1可以直观地看出二元函数在点$\left(1/2,1/2\right)$ 处的偏导数的几何意义，如图中所示切线的斜率$k=-2$ ，${z^{\prime }}_x\left(1/2,1/2\right)=-2$ 。

例2 函数$f\left(x,y\right)=x^3-y^3+3x^2+3y^2-9x$ 的极值问题的可视化。

分析 根据函数极值的必要条件与充分条件即可求出极值点与相应的极值。

根据必要条件，易求四个驻点：$\left(1,0\right),\left(1,2\right),\left(-3,0\right),\left(-3,2\right)$ ；

根据充分条件，验证这四个驻点中极大值点为$\left(-3,2\right)$ ，极小值点为$\left(1,0\right)$ 。

我们也可以作图观察等高线的形状，确定极值点所在的位置。

代码如下：

x=-4:0.1:2;

y=-2:0.1:4;

[X,Y]=meshgrid(x,y);

Z=X.^3-Y.^3+3*X.^2+3*Y.^2-9*X;

mesh(X,Y,Z);

contour(X,Y,Z,40)

hold on

scatter(1,0,'filled','k')

hold on

scatter(-3,2,'filled','b')

输出函数图形见图2与图3。

Figure 2. Function image in Example 2

图2. 例2中函数图像

Figure 3. Function contour plot in Example 2

图3. 例2中函数等高线图

从等高线图3可以看出，梯度方向即为等高线的法线方向且指向函数增加的方向。由此可知，极大值点为$\left(-3,2\right)$ ，极小值点为$\left(1,0\right)$ 附近被等高线环绕；而在其他非极值点周围没有等高线环绕，这为我们解决多元函数极值问题提供了一种图形化方法。

3. Matlab可视化二重积分的计算

3.1. 预备知识

二重积分${\displaystyle {\iint }_{D}f\left(x,y\right)\text{d}x\text{d}y}$ 的几何意义为([6]: p. 130)：设有一立体，它的底是xoy平面上的闭区域D，它的侧面是以D的边界曲线为准线而母线平行于z轴的柱面，它的顶是曲面$z=f\left(x,y\right)\left(f\left(x,y\right)\ge 0\right)$ 且在D上连续，这种曲顶柱体的体积可以表示为$V={\displaystyle {\iint }_{D}f\left(x,y\right)\text{d}x\text{d}y}$ 。

3.2. 二重积分求体积Matlab程序实现

由于实际计算重积分时积分上下限难以确定，计算繁琐，空间体V的几何形状也非常复杂，以下我们通过Matlab可视化作出函数图像帮助学生理解利用二重积分的几何意义求曲顶柱体的体积，加深对重积分的几何意义的理解。

例3 求球体$x^2+y^2+z^2\le 16$ 被圆柱面$x^2+y^2=4x$ 所截得的立体体积V。

分析 将曲面一般式方程改写为参数形式方程：

球面方程$x^2+y^2+z^2=16$ 的参数形式$\{\begin{array}{l}x=4\sin \phi \cos \theta \\ y=4\sin \phi \sin \theta \\ z=4\cos \phi \end{array}$ ，其中$\{\begin{array}{l}\phi \in \left[0,\pi \right]\\ \theta \in \left[0,2\pi \right]\end{array}$ ；

圆柱面$x^2+y^2=4x$ 的参数形式$\{\begin{array}{l}x=2+2\cos t\\ y=2\sin t\\ z=s\end{array}$ ，其中$\{\begin{array}{l}t\in \left[0,2\pi \right]\\ s\in \left(-\infty ,+\infty \right)\end{array}$ 。

Figure 4. The curved top cylinder of the sphere of example 2 cut off by a cylinder

图4. 例2球体被圆柱面所截得的曲顶柱体

立体图形代码如下：

phi=0:0.1:pi;

theta=-pi:0.1:pi;

[phi,theta]=meshgrid(phi,theta);

x1=4.*sin(phi).*cos(theta);

y1=4.*sin(phi).*sin(theta);

z1=4.*cos(phi);

surf(x1,y1,z1) %绘制球面

hold on

s=-5:0.1:5;

t=-pi:0.1:pi;

[s,t]=meshgrid(s,t);

x2=2+2*cos(t);

y2=2*sin(t);

z2=s;

surf(x2,y2,z2) %绘制圆柱面

球体被圆柱面所截之后形成的几何体如图4，其体积可以转化为第一卦限内以xoy平面$y=\sqrt{4x-x^2}$ 与x轴围成的闭区域作为底面$D_1$ ，以曲面$z=\sqrt{16-x^2-y^2}$ 为顶构成的曲顶柱体体积的4倍，即所求体积$V=4{\displaystyle {\iint }_{D_1}\sqrt{16-x^2-y^2}\text{d}x\text{d}y}$ ，考虑极坐标代换：$x=\rho \cos \theta ,y=\rho \sin \theta$ ，则

$V=4{\displaystyle {\iint }_{D_1}\sqrt{16-x^2-y^2}\text{d}x\text{d}y}=4{\displaystyle {\int }_0^{\frac{\pi }{2}}\text{d}\theta }{\displaystyle {\int }_0^{4\cos \theta }\sqrt{16-{\rho }^2}\rho \text{d}\rho }=\frac{128\pi }{3}-\frac{512}{9}$ .

Figure 5. Figure D₁ of the bottom surface of the curved top column truncated by example 2

图5. 例2所截得曲顶柱体底面D₁图形

平面区域底面$D_1$ 图形(见图5)代码如下：

x=0:0.1:4;

y=sqrt(4*x-x.^2);

plot(x,y);

hold on

area(x,y);

三重积分计算代码如下：

syms rho theta

f=sqrt(16-rho^2)*rho;

I=4*int(int(f,rho,0,4*cos(theta)),theta,0,pi/2)

复制以上代码至“实时编辑器”窗口，运行即可。

4. Matlab可视化三重积分的计算

4.1. 预备知识

三重积分${\displaystyle {\iiint }_{\Omega }f\left(x,y,z\right)\text{d}x\text{d}y\text{d}z}$ 在空间中已不具有明显的几何意义，我们可以从物理意义上理解为空间体积V上物体质量([6]: p. 154)，即：$f\left(x,y,z\right)$ 表示某物体$\left(x,y,z\right)$ 处的密度，Ω是该物体所占有的空间闭区域，$f\left(x,y,z\right)$ 在Ω上连续，那么该物体的质量为$m={\displaystyle {\iiint }_{\Omega }f\left(x,y,z\right)\text{d}x\text{d}y\text{d}z}$ 。

4.2. 三重积分Matlab程序实现

通过例题，利用Matlab作出函数图像帮助学生理解三重积分的计算。

例4 计算三重积分${\displaystyle {\iiint }_{\Omega }\sqrt{x^2+y^2}\text{d}x\text{d}y\text{d}z}$ ,其中Ω是由一抛物面$z_1=18-\frac{1}{2}\left(x^2+y^2\right)$ 及其另一抛物面$z_2=5+x^2+y^2$ 所包围的闭区域。

Figure 6. The closed region surrounded by two paraboloids in Example 4

图6. 例4中两抛物面所包围的闭区域

分析 两抛物面所围成的闭区域可以通过Matlab作出，代码如下：

z1 = @(x,y)x.^2*(-0.5)-y.^2*0.5+18;

z2 = @(x,y)x.^2+y.^2+5;

z3 = @(x,y)x.^2*0+y.^2*0;

fun = @(x,y)z1(x,y)-z2(x,y);

x = -3:.25:3; y = -3:.25:3;

[X,Y] = meshgrid(x,y);

Z1 = z1(X,Y);

surf(X,Y,Z1); %绘制z1抛物面

hold on;

Z2 = z2(X,Y);

surf(X,Y,Z2); %绘制z2抛物面

hold on;

h = fimplicit(fun);

h.MeshDensity = 25;

x2 = h.XData; y2 = h.YData;

Z = z2(x2,y2);

plot3(x2,y2,Z,'r*'); %绘制两抛物面交线

Z3 = z3(x2,y2);

fill3(x2,y2,Z3,'c');

for i = 1:length(x2)

a = x2(i)^2+y2(i)^2+5.;

plot3([x2(i),x2(i)],[y2(i),y2(i)],[0,a],'b'); %绘制两抛物面所围立体的投影

hold on;

end

由图6可知，上抛物曲面为$z=18-\frac{1}{2}\left(x^2+y^2\right)$ ，下抛物曲面为$z=5+x^2+y^2$ ，其交线为$\{\begin{array}{l}{x}^2+y^2=9\\ z=14\end{array}$ ，在xoy平面上的投影为$x^2+y^2\le 9$ . 根据柱坐标变换$\{\begin{array}{l}x=\rho \cos \theta \\ y=\rho \sin \theta \\ z=z\end{array}$ ，得

${\displaystyle {\iiint }_{\Omega }\sqrt{x^2+y^2}\text{d}x\text{d}y\text{d}z}={\displaystyle {\int }_0^{2\pi }\text{d}\theta }{\displaystyle {\int }_0^3\text{d}r}{\displaystyle {\int }_{5+r^2}^{18-\frac{1}{2}{r}^2}{r}^2\text{d}z}=\frac{441\pi }{5}$ .

三重积分计算代码如下：

syms r theta z

I=int(int(int(r^2,z,5+r^2,18-0.5*r^2),r,0,3),theta,0,2*pi)

Matlab是一种高性能的数值计算和科学计算软件，它提供了强大的符号计算和数值积分功能。通常的Matlab的符号计算功能就能够计算二重积分和三重积分，如果涉及稍复杂的函数或积分区域时，Matlab中数值积分函数integral2与integral3对于小规模的积分计算问题也可以解决。下面再举几例利用Matlab计算重积分，展示其强大的计算功能。

例5 计算${\displaystyle {\iiint }_{\Omega }{z}^2}\text{d}x\text{d}y\text{d}z,\Omega :\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1$ 。

分析：需将三重积分化为三个定积分，再通过代码实现计算：

Step 1 考虑在此椭球面上，用$z=h$ 截痕，则z的范围为$\left[-c,c\right]$ ；

Step 2 又在截面上，椭圆为$D:\frac{x^2}{a^2}+\frac{y^2}{b^2}=1-\frac{z^2}{c^2}$ ；

Step 3 该三重积分可以化为：${\displaystyle {\int }_{-c}^c\left({\displaystyle {\iint }_D{z}^2\text{d}x\text{d}y}\right)\text{d}z}={\displaystyle {\int }_{-c}^c\left(\pi ab\left(1-\frac{z^2}{c^2}\right)z^2\right)\text{d}z}$ 。

最后一步用代码实现计算结果即可，代码如下：

syms a b c positive

syms theta r x y z

x=a*sqrt(1-z^2/c^2)*r*cos(theta);

y=b*sqrt(1-z^2/c^2)*r*sin(theta); %广义极坐标变换

J=jacobian([x y],[r,theta]); %雅可比矩阵

jacobi_j=det(J); %雅可比行列式

p1=int(jacobi_j,r,0,1);

p2=int(p1,theta,0,2*pi);

p3=int(p2*z^2,z,-c,c)

输出结果为$\frac{4}{15}\pi ab{c}^3$ 。

例6 计算${\displaystyle {\iiint }_{\Omega }\frac{{\text{e}}^z}{\sqrt{x^2+y^2}}\text{d}x\text{d}y\text{d}z},\Omega :z=\sqrt{x^2+y^2},z=1,z=2$ .

解 Step 1 考虑在此空间锥面上，截痕$z=h$ ，则z的范围为$\left[1,2\right]$ ；

Step 2 在截痕面上，平面区域圆为$D:x^2+y^2=z^2$ ；

Step 3 该三重积分可化为：${\displaystyle {\int }_1^2\left({\displaystyle \underset{D}{\iint }\frac{{\text{e}}^z}{\sqrt{x^2+y^2}}\text{d}x\text{d}y}\right)\text{d}z}={\displaystyle {\int }_1^2\text{d}z\left({\displaystyle \underset{D}{\iint }\frac{1}{\sqrt{x^2+y^2}}\text{d}x\text{d}y}\right)}$ ，

最后一步用代码实现计算结果即可，代码如下：

syms x y z theta r

x=r*cos(theta);

y=r*sin(theta);

J=jacobian([x,y],[r,theta]);

jacobi_j=det(J);

f=1/sqrt(x^2+y^2)*jacobi_j;

p1=int(f,r,0,z);

p2=int(p1,theta,0,2*pi);

p3=int(exp(z)*p2,z,1,2)

输出结果为$2\pi {\text{e}}^2$ 。

在运用多元积分解决具体实际问题时，首先要通过分析问题，列出具体积分式，然后将多元积分转化为几个定积分的计算形式，最后利用Matlab程序代码求解。

5. 结论

本文所给出的利用Matlab计算的例题均要求多元函数在指定有界区域内是连续函数。另外，值得注意的是，如计算积分${\displaystyle {\int }_0^1\left({\displaystyle {\int }_x^1{\text{e}}^{y^2}}\text{d}y\right)\text{d}x}$ 时，需要我们先交换积分次序${\displaystyle {\int }_0^1\left({\displaystyle {\int }_x^1{\text{e}}^{y^2}}\text{d}y\right)\text{d}x}={\displaystyle {\int }_0^1\left({\displaystyle {\int }_0^y{\text{e}}^{y^2}}\text{d}x\right)\text{d}y}$ ，然后再计算才可以得到结果$\frac{1}{2}\left(\text{e}-1\right)$ ；但对于Matlab程序代码而言，无需考虑是否先要交换顺序，Matlab具有自适应计算性，都可以得到正确的结果。

Matlab在微积分绘图与计算方面的功能非常强大，在数学类课程教学过程中，引入软件实践增强图形的直观性，计算的快捷性；有利于培养学生对数学学习的兴趣。特别在多变量微积分学的教学过程中，Matlab软件实践与课堂教学相结合弥补了传统教学中的不足。软件技术已经渗透到各个领域，教学与软件实践相结合更好地提高了教学效果。

基金项目

浙江省普通本科高校“十四五”教学改革项目(jg20220747)；高等学校大学数学教学研究与发展中心2023年教学改革项目(CMC20230301)；浙江理工大学2023年教育教学改革研究项目(jgzd202301)。

NOTES

^*通讯作者。