多元微积分实验.ppt

《多元微积分实验.ppt》由会员分享，可在线阅读，更多相关《多元微积分实验.ppt（59页珍藏版）》请在三一办公上搜索。

1、3.1 曲面绘图,3.2 多元函数微分,3 多元微积分实验,3.3 多元函数积分3.4 常微分方程求解,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=mes

2、hgrid(a,b)将向量a,b生成axb的网格矩阵x,y。,实验目的 学习和掌握空间曲面的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);%s

3、urf(x,y,z);,例3.1.2 画出曲面 在矩阵区域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

4、:n条等高线;h:按向量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);figur

5、e(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);t=0:0.01:pi+0

6、.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 onh1=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;camlight(-220,-170);axis equa

7、l,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:2*pi+0.01;s=t;x=2*sin(t)*cos(s);y=2*sin(t)*sin

8、(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 onh1=surf(x1,y1,z1);,%修饰hold off;view(140,9);light(position,2 1 2);axis equal;shading interp;axis off;camlight(-220,-170);set(h,facecolor,0,0.8,0)

9、;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 onh1=surf(x1,y1,z1);view(140,9);light(position,2 1 2);lighting ph

10、ong;shading interp;%axis off;camlight(-220,-170);axis equal,3.2 多元函数微分,3.2.1 多元函数极限 多元函数沿不同路径趋于某一点时,其极限可能会出现不同的结果.实验目的 学习和掌握用matlab工具求解多元函数极限问题.例3.2.1 求syms x ylimit(limit(x2+y2)/(sin(x)+cos(y),0),pi)limit(limit(1-cos(x2+y2)/(x2+y2)*exp(x2*y2),0),0),3.2.2 多元函数偏导数及全微分 多元函数对某一变量的偏导数可以利用一元函数导数的命令.实验目的

11、学习和掌握用matlab工具求解多元函数的偏导数和全微分.多元函数的偏导数和全微分问题化为一元函数导数及微分问题.,例3.2.2 设 的一阶偏导数和全微分及syms x y dx dyz=atan(x2*y);zx=diff(z,x)zy=diff(z,y)dz=zx*dx+zy*dyzxx=diff(zx,x)zxy=diff(zx,y),例设 f(x,y,z)=xy2+yz2+zx2 求：求fxx(0,0,1),fxz(1,0,2),fyz(0,-1,0)及 fzzx(2,0,1)syms x y zf=x*y2+y*z2+z*x2;fxz=diff(diff(f,x),z);x=1;y=

12、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)给出(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),

13、Y(:,11:21),Z(:,11:21);view(120,18),3.2.3.2 切线与法平面 空间曲线处的切线的方向向量为用matlab命令得到切向量：过点 的切线方程F为：即：,转化为matlab语句为：即：过点 的法平面方程G为：即：转化为matlab语句为：,例3.2.4 空间曲线处的切线方程和法平面方程,并画图表示.syms t x y zx1=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 tF=-x;y;z+x1

14、;y1;z1+s0*tG=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)/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 onplot3(x2,

15、y2,z2);hold onmesh(x,y,z);axis equalview(-45,15),另法:的切线与法平面.切线:法平面:得:,Smys t x y zX1=;y1=;z1=;Xt=diff(x,t);T=t0;x0=;Xt0=subs(xt);Syms t x y z%切线方程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 onplot3

16、(x2,y2,z2);hold onmesh(x,y,z);axis equalview(-45,15),syms t x y zx1=3*sin(t);y1=3*cos(t);z1=5*t;xt=diff(x1,t);yt=diff(y1,t);zt=diff(z1,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 z2x2=x0+xt0*tt,y2=y0+yt0*tt,z2=z0+zt0*ttdisp(法平面方程)z=-(x

17、-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);y2=subs(y2);z2=subs(z2);z=subs(z);plot3(x1,y1,z1);hold onplot3(x2,y2,z2);hold onmesh(x,y,z);axis equalview(-45,15),3.2.3.3 切平面与法线空间曲面处的切平面法向量:过点 切平面方程:法线方程,例3.2.5 设曲面方程处的切平面

18、和法线方程.syms x y z tf=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,z0);disp(切平面方程：);F=x-x0,y-y0,z-z0*ndisp(法线方程:);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=s

19、ubs(fy);fz=subs(fz);syms t x3 y3 z3x0=1;y0=1;z0=4;disp(切平面方程)%x3,y3,z3z3=-(fx*(x3-x0)+fy*(y3-y0)-fz*z0)/fzdisp(法线方程)%x2,y2,z2x2=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.

20、2:2,-2:0.2:2);z3=subs(z3);mesh(x1,y1,z1);hold onplot3(x2,y2,z2);hold onmesh(x3,y3,z3);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 onquiver(x

21、,y,px,py),hold off,3.2.4 多元函数的极值X=fminsearch(ff,x0)%单纯形法X=fminunc(ff,x0)%拟牛顿法实验目的 学习和掌握用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工

22、具求解二重积分问题.3.3.1.1 积分限为常数例3.3.1 计算disp(符号法);syms x ys=vpa(int(int(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 onplot(x,y2,2.5,-0.5:0.01:3);text(1.5,0.2,y=1/2x);text(1.5,2,y=(2x)1

23、/2);axis(-0.5 3-0.5 3);,%求交点syms x yy1=2*x*y-1;y2=y-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 zs=vpa(int(int(int(x*y*z,x,0,1),y,1

24、,2),z,2,3),disp(数值法);w=inline(x*y*z,x,y,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 onx1,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

25、(t);y=(0*s+2)*sin(t);z=s*(0*t+2);mesh(x,y,z);,%确定积分限syms x y zf1=(z=8-x2-y2);f2=(x2+y2=4);x1,y1,z1=solve(f1,f2,x,y,z)%分步积分syms x y zf=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进行符号求解.

26、在matlab中数值法有一些常用解法: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*l

27、og(x)+y=a*x*(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 常微分方程的特解 如果给定微分方程满足初始条件则求方程的特解的调用格式为:dsolv

28、e(eqn,condition1,condition n,var),例3.4.3 求下列常微分方程给定初始条件的特解y=dsolve(D2y=cos(2*x)-y,y(0)=1,Dy(0)=0,x)simplify(y),3.4.2 常微分方程的数值求解 初值问题数值求解的matlab常用格式 t,y=oder45(odefun,tspan,y0)4/5b阶龙格-库塔-芬尔格算法 odefun表示f(t,y)的函数句柄或inline函数,t是标量,y是标量或向量.tspan若为t0,tf,表示自变量初值与终值,若为t0,t1,tn表示输出节点列向量.y0表示初值向量y y:输出数值解矩阵 t输

29、出参数t,实验目的 学习和掌握用matlab工具解决常微分方程数值求解问题例3.4.4 求微分方程的解析解和数值解,并画出对应图形进行比较.clear;解析解y=dsolve(Dy=y-2*x/y,y(0)=1,x)x=0:0.1:4;y=subs(y);fun=inline(s-2*t/s,t,s);t,s=ode45(fun,0,4,1);disp(数值解);t,splot(x,y,o-,t,s,.);,例3.4.5 求解微分方程令:则原微分方程化为一阶微分方程组:,function dy=weifen1(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1

30、-y(1)2)*y(2)-y(1);输入命令:t,y=ode15s(weifen1,0 3000,0 1);plot(t,y(:,1),-);,例3.4.6 求刚性微分方程的解f(1)=inline(-0.01*y(1)-99.99*y(2);f(2)=inline(-100*y(2);t,y=ode15s(weifen2,0 400,2 1);plot(t,y,-);,输入命令:t,y=ode15s(weifen2,0 400,2 1);plot(t,y,-);另法:f(1)=inline(-0.01*y(1)-99.99*y(2);f(2)=inline(-100*y(2);t,y=ode15s(weifen2,0 400,2 1);plot(t,y,-);,