多元微积分实验

1、3 多元微积分实验,3.1 曲面绘图,曲面的一般方程是F(x,y,z)=0，在matlab中将曲面的点(x,y,z)的坐标先表示出来，再使用对应的曲面绘图函数。matlab常用的绘图函数有： plot3,mesh,surf等,1）plot3 一组平行平面上截线（曲线族）方式来表示曲面。 （2）mesh 两组相交的平行平面上网格线方式来表示曲面。 （3）surf 两组网状线和补片填充色彩的方式表示曲面。 （4）meshc 并带有等高线。 （5）surfc 并带有等高线。 （6）surfl 并带有阴影。 （7）x,y=meshgrid(a,b) 将向量a,b生成axb的网格矩阵x,y,实验目的 学

2、习和掌握空间曲面的matlab作图方法。 （1）平面网格坐标矩阵的生成 x=a:dx:b; y=c:dy:d; X,Y=meshgrid(x,y); %利用meshgrid生成 (2)函数 z表达式(点运算) （3）绘图函数绘图 绘图函数（x,y,z） x,y=meshgrid(-9:0.3:9); %ac相同,bd相同,例3.1.1 画出曲面 在矩形区域 D:-4 b=-5:0.2:5; x,y=meshgrid(a,b); z=exp(-x.2-y.2)/2)/2/sqrt(2*pi); plot3(x,y,z); % mesh(x,y,z); % surf(x,y,z,例3.1.2 画出

3、曲面 在矩阵区域D：-9 plot3(x,y,z); % mesh(x,y,z); % surf(x,y,z,当区域 时 clear r=0:0.1:9;t=0:pi/50:2*pi; R,T=meshgrid(r,t); x=R.*cos(T);y=R.*sin(T); z=sin(sqrt(x.2+y.2)./sqrt(x.2+y.2+eps); plot3(x,y,z); % mesh(x,y,z); % surf(x,y,z,3.1.2 等高线的绘制 在matlab中绘制等高线常用的函数: contour(x,y,z,选项) 画出在xoy平面上的曲面的等高线. 选项n:n条等高线;h:

4、按向量h的值画出投影在xoy平面上的曲面的等高线;str:所指定颜色与线形. contourf(x,y,z) 填充颜色 c,h=contour(z) z曲面上的点,c为xoy平面上的等高线阵,h高度列向量,Clabel(c,h) Clabel(c,manual) Contour3(x,y,x,n) 在空间中画,例3.1.3 画出曲面 的各种等高线. clear;clf x,y=meshgrid(-4:0.2:4); z=x.3+y.3-12*x-12*y; figure(1);mesh(x,y,z); figure(2); c,h=contour(x,y,z); clabel(c,h); fi

5、gure(3); h1=-28 -16 -8 0 6 18 26; c1=contour(z,h1); clabel(c1); figure(4);contourf(z); contour3(z,10,例3.1.4 画出两圆柱面 的相交图形. figure(1); t=0:0.03:2*pi;s=-2:0.03:2; x=(0*s+1)*cos(t);y=(0*s+1)*sin(t);z=s*(0*t+1); mesh(x,y,z); hold on; mesh(x,z,y); axis equal,(0*s+1)产生s向量大小的向量,figure(2); figure(color,1,1,1

6、); t=0:0.01:pi+0.01;s=0:0.01:2; x=cos(t)*(0*s+1); y=sin(t)*(0*s+1); z=(0*t+1)*s; h=surf(x,y,z); hold on h1=surf(x,-y,-z);h2=surf(x,-y,z); h3=surf(x,y,-z);h4=surf(z,y,x); h5=surf(z,-y,x);h6=surf(-z,y,x); h7=surf(-z,-y,x); view(-128,23); light(position,2 1 2); lighting phong; shading interp; %axis off

7、; camlight(-220,-170); axis equal,figure(1); t=0:0.03:2*pi;s=-2:0.03:2; x=(0*s+1)*cos(t);y=(0*s+1)*sin(t);z=s*(0*t+1); mesh(x,y,z); hold on; mesh(x,z,y); view(-128,23); light(position,2 1 2); lighting phong; shading interp; %axis off; camlight(-220,-170); axis equal,例3.1.5 画出圆柱面 和球面 的相交图形. t=0:0.01:

8、2*pi+0.01;s=t; x=2*sin(t)*cos(s);y=2*sin(t)*sin(s);z=2*cos(t)*(0*s+1); t1=0:0.01:2*pi+0.01;s1=-2:0.01:+2; x1=(1+cos(t1)*(0*s1+1);y1=sin(t1)*(0*s1+1); z1=(0*t1+1)*s1; figure(color,1,1,1); h=surf(x,y,z); hold on h1=surf(x1,y1,z1,修饰 hold off;view(140,9); light(position,2 1 2);axis equal; shading interp

9、;axis off;camlight(-220,-170); set(h,facecolor,0,0.8,0); set(h1,facecolor,1,0,1,t=0:0.01:2*pi+0.01;s=t; x=2*sin(t)*cos(s);y=2*sin(t)*sin(s);z=2*cos(t)*(0*s+1); t1=0:0.01:2*pi+0.01;s1=-2:0.01:+2; x1=(1+cos(t1)*(0*s1+1);y1=sin(t1)*(0*s1+1); z1=(0*t1+1)*s1; figure(color,1,1,1); h=surf(x,y,z);hold on h1

10、=surf(x1,y1,z1); view(140,9);light(position,2 1 2); lighting phong; shading interp;%axis off; camlight(-220,-170); axis equal,3.2 多元函数微分,3.2.1 多元函数极限 多元函数沿不同路径趋于某一点时,其极限可能会出现不同的结果. 实验目的 学习和掌握用matlab工具求解多元函数极限问题. 例3.2.1 求 syms x y limit(limit(x2+y2)/(sin(x)+cos(y),0),pi) limit(limit(1-cos(x2+y2)/(x2+

11、y2)*exp(x2*y2),0),0,3.2.2 多元函数偏导数及全微分 多元函数对某一变量的偏导数可以利用一元函数导数的命令. 实验目的 学习和掌握用matlab工具求解多元函数的偏导数和全微分. 多元函数的偏导数和全微分问题化为一元函数导数及微分问题,例3.2.2 设 的一阶偏导数和全微分及 syms x y dx dy z=atan(x2*y); zx=diff(z,x) zy=diff(z,y) dz=zx*dx+zy*dy zxx=diff(zx,x) zxy=diff(zx,y,例设 f (x, y, z) = xy2 + yz2 + zx2 求 ： 求fxx( 0, 0, 1)

12、, fxz (1, 0, 2 ), fyz ( 0, -1, 0)及 fzzx ( 2, 0, 1) syms x y z f=x*y2+y*z2+z*x2; fxz=diff(diff(f,x),z); x=1;y=0;z=2; fxz=subs(fxz) %求函数值,3.2.3 微分法在几何上的应用 这里主要讨论微分法在法线、切线与法平面、切平面与法线、数值梯度上的应用。 实验目的 学习和掌握用matlab工具来解决微分法在几何上的应用。 3.2.3.1 法线 surfnorm(X,Y,Z) 绘制(X,Y,Z) 所表示的曲面的法线。 Nx,Ny,Nz=surfnorm(X,Y,Z) 给出(

13、X,Y,Z) 所表示的曲面的法线数据,例3.2.3 绘制半个椭圆 绕x轴旋转一周得到的曲面的法线。 y=-1:0.1:1;x=2*cos(asin(y); X,Y,Z=cylinder(x,20); surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21); view(120,18,3.2.3.2 切线与法平面 空间曲线 处的切线的方向向量为 用matlab命令得到切向量： 过点 的切线方程F为： 即,转化为matlab语句为： 即： 过点 的法平面方程G为： 即： 转化为matlab语句为,例3.2.4 空间曲线 处的切线方程和法平面方程,并画图表示. syms t

14、 x y z x1=3*sin(t);y1=3*cos(t);z1=5*t; s1=jacobian(x1,y1,z1,t); t=pi/4; x0=3*sin(t);y0=3*cos(t);z0=5*t; s0=subs(s1); syms t F=-x;y;z+x1;y1;z1+s0*t G=x-x0,y-y0,z-z0*s0,切线方程: 法平面,t=-pi:0.1:pi;x,y=meshgrid(-3:0.2:3);tt=-3:0.2:3; x1=3*sin(t);y1=3*cos(t);z1=5*t; x2=3*sqrt(2)/2+3*sqrt(2)*tt/2; y2=3*sqrt(2

15、)/2-3*sqrt(2)*tt/2;z2=5*pi/4+5*tt; z=(3*sqrt(2)/2*(x-3*sqrt(2)/2)-3*sqrt(2)/2*(y-3*sqrt(2)/2)-25*pi/4)/(-5); plot3(x1,y1,z1);hold on plot3(x2,y2,z2);hold on mesh(x,y,z); axis equal view(-45,15,另法: 的切线与法平面. 切线: 法平面: 得,Smys t x y z X1=;y1=;z1=; Xt=diff(x,t); T=t0; x0=; Xt0=subs(xt); Syms t x y z %切线方程

16、 X1=x0+xt0*t,y1=y0+xt0*t,z1=z0+zt0*t %法平面方程 Z=-(x-x0)*xt0+(y-y0)*yt0-z0)/zt0,画图 T=;x,y=meshgrid();tt= X1=subs(x1);y1=y1;z1=z1; Z=z; plot3(x1,y1,z1);hold on plot3(x2,y2,z2);hold on mesh(x,y,z); axis equal view(-45,15,syms t x y z x1=3*sin(t);y1=3*cos(t);z1=5*t; xt=diff(x1,t);yt=diff(y1,t);zt=diff(z1,

17、t); t=pi/4; x0=subs(x1);y0=subs(y1);z0=subs(z1); xt0=subs(xt);yt0=subs(yt);zt0=subs(zt); disp(切线方程) syms tt x2 y2 z2 x2=x0+xt0*tt,y2=y0+yt0*tt,z2=z0+zt0*tt disp(法平面方程) z=-(x-x0)*xt0+(y-y0)*yt0-z0)/zt0,画图 t=-pi:0.1:pi;x,y=meshgrid(-3:0.2:3);tt=-3:0.2:3; x1=subs(x1);y1=subs(y1);z1=subs(z1); x2=subs(x2

18、);y2=subs(y2);z2=subs(z2); z=subs(z); plot3(x1,y1,z1);hold on plot3(x2,y2,z2);hold on mesh(x,y,z); axis equal view(-45,15,3.2.3.3 切平面与法线 空间曲面 处的切平面法向量: 过点 切平面方程: 法线方程,例3.2.5 设曲面方程 处的切平面和法线方程. syms x y z t f=3*x2+y2-z;w=x,y,z; s=jacobian(f,w);x0=1;y0=1;z0=4; n=subs(s,x,x0);n=subs(n,y,y0);n=subs(n,z,z

19、0); disp(切平面方程：);F=x-x0,y-y0,z-z0*n disp(法线方程:); G=-x;y;z+x0,y0,z0+n*t %G=-x,y,z+x0,y0,z0+n*t,syms t x y z %x,y,z曲面 f=3*x2+y2-z; fx=diff(f,x);fy=diff(f,y);fz=diff(f,z); x=1;y=1;z=4; fx=subs(fx);fy=subs(fy);fz=subs(fz); syms t x3 y3 z3 x0=1;y0=1;z0=4; disp(切平面方程) %x3,y3,z3 z3=-(fx*(x3-x0)+fy*(y3-y0)-

20、fz*z0)/fz disp(法线方程) %x2,y2,z2 x2=x0+fx*t,y2=y0+fy*t,z2=z0+fz*t,u=0:0.1:1.5;v=0:0.1:2*pi;t=-1:0.1:0.5; x1=u*cos(v);y1=sqrt(3)*u*sin(v);z1=3*u.2*(0*v+1); x2=subs(x2);y2=subs(y2);z2=subs(z2); x3,y3=meshgrid(0:0.2:2,-2:0.2:2); z3=subs(z3); mesh(x1,y1,z1);hold on plot3(x2,y2,z2);hold on mesh(x3,y3,z3);

21、view(156,68);axis equal,3.2.3.4 数值梯度 FX,FY=gradient(F,h) 求二元函数的梯度,FX,FY为沿x,y方向的数值导数. 例3.2.6 已知二元函数 步长为0.2,试求梯度向量并画图,v=-2:0.2:2; x,y=meshgrid(v); z=x.*exp(-x.2-y.2); px,py=gradient(z,2,0.2); contour(x,y,z);hold on quiver(x,y,px,py),hold off,3.2.4 多元函数的极值 X=fminsearch(ff,x0) %单纯形法 X=fminunc(ff,x0) %拟牛

22、顿法 实验目的 学习和掌握用matlab工具求解多元函数的极值问题. 例3.2.7 求二元函数 在原点附近的极大值 (问题等价于一元函数的极小值) fun=inline(x(1)4+x(2)4-4*x(1)*x(2)-5); x,y=fminsearch(fun,0,0,3.3 多元函数积分,3.3.1 二重积分 二重积分可以转化为二次积分运算, 即 或 可以用matlab来计算两个定积分, 还可用dblquad进行数值计算,实验目的 学习和掌握用matlab工具求解二重积分问题. 3.3.1.1 积分限为常数 例3.3.1 计算 disp(符号法); syms x y s=vpa(int(i

23、nt(xy,x,0,1),y,1,2) disp(数值法); fun=inline(x.y,x,y); s=dblquad(fun,0,1,1,2,3.3.1.2 内积分限为函数 例3.3.2 计算 ,其中Dxy由曲线 所围成的平面区域. x=0.01:0.01:3; y1=1./(2*x);y2=sqrt(2*x); plot(x,y1);hold on plot(x,y2,2.5,-0.5:0.01:3); text(1.5,0.2,y=1/2x); text(1.5,2,y=(2x)1/2); axis(-0.5 3 -0.5 3,求交点 syms x y y1=2*x*y-1;y2=y

24、-sqrt(2*x); x0,y0=solve(y1,y2); %求积分 f=exp(-x2-y2); y1=1/(2*x);y2=sqrt(2*x); jfy=int(f,y,y1,y2); jfx=int(jfy,x,x0,2.5); jf2=vpa(jfx,3.3.2 三重积分 三重积分可化为三次积分计算: 实验目的 学习和掌握用matlab来计算三重积分问题. 3.3.2.1 积分限为常数 例3.3.3 计算 %符号法 syms x y z s=vpa(int(int(int(x*y*z,x,0,1),y,1,2),z,2,3,disp(数值法); w=inline(x*y*z,x,y

25、,z); s=triplequad(w,0,1,1,2,2,3) 3.3.2.2 内积分限为函数 例3.3.4 计算 ,其中积分区域V是由旋转抛物面 所围成的空间闭区域,画区域图 t,r=meshgrid(0:0.05:2*pi,0:0.05:2); x=r.*cos(t);y=r.*sin(t);z=8-x.2-y.2; mesh(x,y,z);hold on x1,y1,z1=cylinder(2,30); z2=4*z1;mesh(x1,y1,z2) t=0:0.03:2*pi;s=-2:0.03:2; x=(0*s+2)*cos(t);y=(0*s+2)*sin(t);z=s*(0*t

26、+2); mesh(x,y,z,确定积分限 syms x y z f1=(z=8-x2-y2);f2=(x2+y2=4); x1,y1,z1=solve(f1,f2,x,y,z) %分步积分 syms x y z f=x+exp(y)+sin(z);z1=0;z2=8-x2-y2; jfz=int(f,z,z1,z2); jfx=int(jfz,x,x1(1),x1(2); jfy=int(jfx,y,-2,2); jf3=vpa(jfy,3.4 常微分方程求解,常微分方程有解析解法和数值解法两种. 在matlab中解析法可用dsolve进行符号求解. 在matlab中数值法有一些常用解法:

27、Euler,Runger-Kutta等,3.4.1 常微分方程(组)符号求解 实验目的 学习和掌握用matlab工具解决常微分方程符号求解问题. 3.4.1.1 常微分方程(组)的通解 利用dsolve命令可以很方便的求得常微分方程 的通解和满足给定条件的特解。但须注意在建立方程时应分别输入为Dy,D2y,D3y,且一般需要指明自变量,主要调用格式: dsolve(eqn,var) %eqn为常微分方程,var为变量,默认t dsolve(eqn1,eqn2,var1,var2,) %常微分方程组 例3.4.1 求下列常微分方程的通解,y1=dsolve(x*Dy*log(x)+y=a*x*(

28、log(x)+1),x) y2=dsolve(D2y+2*Dy+5*y=sin(2*x),x) y3=dsolve(y4-3*x2)*Dy+x*y=0,x) %y3=dsolve(x*y*Dy+x4-3*y2,x) y4=dsolve(D3y+D2y-2*Dy=x*(exp(x)+4),x) simplify(y4,例3.4.2 求下列常微分方程的通解 %x,y=dsolve(Dx=y,Dy=-x) s=dsolve(Dx=y,Dy=-x) y=s.y,x=s.x %s.y 成员运算,3.4.1.2 常微分方程的特解 如果给定微分方程 满足初始条件 则求方程的特解的调用格式为: dsolve(eqn,condition1,