1、1实验五 曲面绘图【实验目的】1 了解二元函数图形的制作。2 空间曲面等高线的制作。3 学习掌握 MATLAB 软件有关的命令。【实验内容】画出函数 的图形, 并画出其等高线。2yxz【实验准备】1曲线绘图的 MATLAB命令MATLAB 中主要用 mesh,surf 命令绘制二元函数图形。mesh(x,y,z) 画网格曲面，这里 x,y,z 是三个数据矩阵，分别表示数据点的横坐标，纵坐标和函数值，该命令将数据点在空间中描出,并连成网格。surf(x,y,z) 画完整曲面, 这里 x,y,z 是三个数据矩阵，分别表示数据点的横坐标，纵坐标和函数值，该命令将数据点所表示曲面画出。可以用 help
2、 mesh, help surf 查阅有关这些命令的详细信息【实验方法与步骤】练习 1 画出函数 的图形,不妨将区域限制在2yxz。用 MATLAB 作图的程序代码为：3,),(yxclear; 2x=-3:0.1:3; %x的范围为-3,3y=-3:0.1:3; %y的范围为-3,3X,Y=meshgrid(x,y); %将向量x,y指定的区域转化为矩阵X,YZ=sqrt(X.2+Y.2); %产生函数值Zmesh(X,Y,Z)结果如图 5.1。图 5.1 是网格线图,如果要画完整的曲面图,只需将上述的MATLAB 代码 mesh(X,Y,Z)改为 surf(X,Y,Z), 结果如图 5.2
3、图 5.1 锥面 图 5.2 锥面要画等高线,需用 contour,contour3 命令. 其中 contour 为二维等高线, contour3 为三维等高线,如画图 5.1 的三维等高线 , MATLAB 代码为：clear; x=-3:0.1:3;y=-3:0.1:3;X,Y=meshgrid(x,y);Z=sqrt(X.2+Y.2);3contour3(X,Y,Z,10) %画10 条等高线xlabel(X-axis),ylabel(Y-axis),zlabel(Z-axis) %三个坐标轴的标记title(Contour3 of Surface) %标题grid on %画网格线结
4、果如图 5.3.-2 02-20201234X-axisY-axisZ-axisContour3 of Surface图 5.3 等高线如画图 5.1 的二维等高线, MATLAB 代码为：clear; x=-3:0.1:3; y=-3:0.1:3;X,Y=meshgrid(x,y); Z=sqrt(X.2+Y.2); contour(X,Y,Z,10)xlabel(X-axis),ylabel(Y-axis)title(Contour of Surface)grid on结果如图 5.4.4-3 -2 -1 0 1 2 3-3-2-10123X-axisY-axisContour of Su
5、rface图 5.4 等高线如果要画 的等高线，则用命令1zclear; x=-3:0.1:3; y=-3:0.1:3;X,Y=meshgrid(x,y); Z=sqrt(X.2+Y.2); contour(X,Y,Z,1 1)结果如图 5.5。-2 0 2-202图 5.5 等高线练习 1 中，函数值 可简单算出。在有些情况下，函数值2yxz不能简单算出。这是因为 x 和 y 的值可能是非均匀间隔的甚至是随),(yxfz5机分布的，也可能使用了不同的坐标系，比如非长方形的网。出现这些情况时,MATLAB 中的函数 griddata 就用来产生经查值后的均匀间隔数据以作图。练习 二次曲面的方程
6、如下 dczbyax22讨论参数 对其形状的影响。cba,本练习的关键在于如何作出三维曲面图形，特别注意在给定 值求yx,时，若有开方运算，一是会出现虚数，二是对实数也有正负两个解。为了使z虚数不出现在绘图中，采用了一种技巧，就是将虚数都换成非数(NaN). MATLAB 代码为：a=input(a=); b=input(b=); c=input(c=);d=input(d=); N=input(N=); %输入参数，N为网格线数目xgrid=linspace(-abs(a), abs(a),N); %建立x网格坐标ygrid=linspace(-abs(b), abs(b),N); %建立y
7、网格坐标x,y=meshgrid(xgrid,ygrid); %确定 个点的x,y网格坐标Nz=c*sqrt(d-y.*y/b2-x.*x/a2); u=1; %u=1，表示z 要取正值z1=real(z); %取z的实部z1for k=2:N-1 %以下7行程序的作用是取消z中含虚数的点for j=2:N-1if imag(z(k,j)=0 z1(k,j)=0; endif all(imag(z(k-1:k+1,j-1:j+1)=0 za(k,j)=NaN; endendend6surf(x,y,z1), hold on %画空间曲面if u=1 z2=-z1; surf(x,y,z2);
8、%u=1时加画负半面axis(-abs(a),abs(a), -abs(b), abs(b), -abs(c), abs(c);endxlabel(x),ylabel(y),zlabel(z)hold off运行程序，当 时的结果见图 5.6, 20,13,45Ndcba当 时的结果见图 5.7,5,i当 时的结果见图 5.8,.0c图 5.6 椭球面7图 5.7 双曲面-505-4-2024-3-2-10123xyz图 5.8 椭球双曲面练习 3 列出求空间两任意曲面的交线的程序。两空间曲面方程连立起来，就形成一个空间曲线的方程。这个曲线能满足两个曲面的方程，因而也就是这两个空间曲面的交线。
9、显示这两个曲面并不难，用两次 mesh 语句即可,但要显示其交线，必须先找到各个交点，因为数值计算得到的是离散点，难以找到两个曲面上完全重合的点，本程序采用了设置限的方法，只要在同一网格点处，两曲面的 z 之之差小于设定限，就认为它是交点，8限值设定几次要才能定的好。下面 MATLAB 程序给出两个空间曲面的交线（当然是空间曲线） ，给出不同的 z1,z2 方程可绘出不同的空间曲线和其交线。x,y=meshgrid(-2:0.1:2,-2:0.1:2); %设定计算和绘图的定义域网格z1=x.2-2*y.2; %第一个曲面方程z2=2*x-3*y; %第二个曲面方程mesh(x,y,z1);
10、hold; mesh(x,y,z2); %再一个图上同时画出两个曲面r0=(abs(z1-z2)zz=r0.*z1; yy=r0.*y; xx=r0.*x; %求这些网格上的坐标值，即交线坐标plot3(xx(r0=0),yy(r0=0),yy(r0=0),*); %画出这些点colormap(gray), hold off %不用彩色而用灰度表示曲面执行此程序得出的曲面见图 5.9.图 5.9 两曲面的交线如果想改表曲面方程，可以在程序中改动第二行和第三行。但这样的程序还不是通用的，最好程序运行时能向用户提问，允许用户输入曲面方程。此时9就要用到字符串功能和 eval 命令。s1=input
11、(输入第一个方程,s);在原来的 z1 方程语句处改为 z1=eval(s1);类似地输入第二个方程。此外，应使用户能给出定义域和间隔。这实现起来比较简单，只要把第一句改为x,y=meshgrid(xmin:dx:xmax,ymin:dy:ymax);其中，xmin，dx ，xmax，ymin，dy ，ymax 可由程序给出屏幕提问，让用户用键盘输入。当然，这样又增加了运行时的麻烦，所以编程时要找一个折衷的选择，要有一定的灵活性又不能太麻烦，应恰到好处。练习 4 用平行界面法讨论由方程 构成的马鞍面形状。2yxz我们只需对练习 3 种的程序作如下修改：定义域网格改为x,y=meshgrid(-
12、10:0.2:10, -10:0.2:10);第一个曲面方程改为z1=(x.2-2*y.2)+eps;第二个曲面（平面）方程改为与z轴正交的水平面，z2=a;为了画z2的曲面图，应使得z2与x,y 有同样的维数，故写为z2=a*ones(size(x);a 可由用户输入,另外用 subplot 把曲面和交线分别画在两张图上 ,并注意把两个分图取成同样比例,便于比较. 因为 z 的范围增大 ,必须把两曲面交点处 z1 和z2 的容差放大到 1.x,y=meshgrid(-10:0.2:10, -10:0.2:10); %设定计算和绘图的定义域网格z1=(x.2-2*y.2)+eps; %第一个曲
13、面方程a=input(a=(-50subplot(1,2,1),mesh(x,y,z1);hold on;mesh(x,y,z2); %分别划出两个曲面v=-10,10,-10,10,-100,100; axis(v), grid %确定第一个分图的坐标系colormap(gray), hold off, %取消彩色,改为灰度r0=abs(z1-z2)zz=r0.*z2; yy=r0.*y; xx=r0.*x; %求这些网格上的坐标值，即交线坐标subplot(1,2,2),plot3(xx(r0=0),yy(r0=0),zz(r0=0),x);%画出交线axis(v), grid %使得第二
14、个分图取第一个分图的坐标系执行此程序,并输入a=8, 得到的三维图形及交线见图 5.10, 当a=-20,得到的三维图形及交线见图5.11,可见从上而下, 其横切面交线发生了很大的变化.图 5.10 马鞍面的水平截面(a=8)11图 5.11 马鞍面的水平截面 (a=-20)练习 5 已经知道曲面上一些点的数据(2,2,80), (3,2,82), (4,2,84), (0,3,79), (2,3,61), (3,3,65), (0,4,84), (1,4,84), (4,4,86), 将这些数据用二元函数插值的方法画出完整的曲面。首先看这些原始数据的柄图，相应的 MATLAB 程序代码为：c
15、lear; x=2,3,4,0,2,3,0,1,4;y=2,2,2,3,3,3,4,4,4;z=80,82,84,79,61,65,84,84,86;stem3(x,y,z); %画火柴杆图命令title(Raw data);xlabel(x),ylabel(y),zlabel(z)结果如图 5.12.12024234050100xRaw datayz图 5.12 柄图显然上面数据是残缺不全的，下面用插值的方法画出完整的曲面，相应的MATLAB 程序代码为：xi=0:0.2:3; yi=2:0.2:4; %选定x,y的范围X,Y=meshgrid(xi,yi); %产生网格向量X,YZ=gri
16、ddata(x,y,z,X,Y,cubic); %cubic采用三角形三次插值mesh(X,Y,Z); title(Griddata);xlabel(x),ylabel(y),zlabel(z)结果如图 5.13.0 12 323460708090xGriddatayz图 5.13 插值曲面13练习 6 （海底测量）表 5-1 给出水面直角坐标(x,y)处水深 z,这时在低潮时测得的。如果船的吃水深度为 5 米，试问在矩形域 中150,2075yx船应避免进入那些区域？表 5-1 水深数据x(m)y(m)z(m)12974140141810828688147818522619513781058
17、58x(m)y(m)z(m)157-69107-8197738145458162-669162844117-389我们首先看测量点的位置：clear; close;x=129 140 108 88 185 195 105 157 107 77 145 162 162 117;y=7 141 28 147 22 137 85 -6 -81 3 45 -66 84 -38;plot(x,y,o);结果如图 5.8.50 100 150 200-100-5005010015014图 5.14 测量点的位置由图 5.8 可见，这是一批不规则数据。由于没有先验函数，我们使用插值法。为了使结果更直观，考虑
18、将 z 的数据转化为相对于海面的高度。相应的MATLAB 程序代码为：z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;h=-z; %数据转化为相对于海面的高度xi=75:5:200; yi=-50:10:150;X,Y=meshgrid(xi,yi);H=griddata(x,y,h,X,Y,cubic);mesh(X,Y,H); view(-60,30); %改变视点结果如图 5.15图 5.15 海底地形图由图 5.15 可见，在(129,7.5)和(162,84)附近各有一块暗礁。进一步，求水深不到 5 米的两个危险区域：15contour(X,Y,H, -5,-5,k) %
19、k表示等高线的颜色为黑色80 100 120 140 160 180 200-50050100150图 5.16 两个危险区域 【练习与思考】1. 画出空间曲线 在 范围内的图形，并画21sin0yxz30,yx出相应的等高线。x,y=meshgrid(-30:0.1:30);z=10*sin(sqrt(x.2+y.2)./sqrt(1+x.2+y.2);meshc(x,y,z)2. 根据给定的参数方程，绘制下列曲面的图形。（1） 椭球面 uzvuyvux sin,cos2,sinco3 x,y,z=ellipsoid(0,0,0,3,2,1,100);16 mesh(x,y,z)axis e
20、qual（2） 椭圆抛物面 24,cos2,sin3uzvyvux syms u v ezmesh(3*u*sin(v),2*u*cos(v),4*u2)（3） 单叶双曲面 uzvuyvux tan4,cose2,sinec3 t,z=meshgrid(-2*pi:pi/10:2*pi,-10:.5:10); x=3*sqrt(z/4).2+1)*sin(t); y=2*sqrt(z/4).2+1)*cos(t); mesh(x,y,z)17（4） 双曲抛物面 3,2vuzyx x,y=meshgrid(-10:0.1:10); z=(x.2-y.2)/3; mesh(x,y,z)（5） 旋转
21、面 uzvyvux,cosln,silnt,z=meshgrid(-2*pi:0.1:2*pi,0.1:0.1:10);x=log(z).*sin(t); y=log(z).*cos(t); mesh(x,y,z)18（6） 圆锥面 uzvyux,cos,sin z,t=meshgrid(-10:0.1:10,-2*pi:0.1:pi); x=z.*sin(t); y=z.*cos(t); mesh(x,y,z)（7） 环面 vzvuyvux sin4.0,si)co4.03(,cos)4.03( syms u v x=(3+0.4*cos(u)*sin(v); y=(3+0.4*cos(u)
22、*cos(v); z=sin(u)*3; z=0.4*sin(u); ezmesh(x,y,z)19（8） 正螺面 vzuyvx4,cos,sin syms u v x=u*sin(v); y=u*cos(v); z=4*u; ezmesh(x,y,z)3. 在一丘陵地带测量搞程,x 和 y 方向每隔 100 米册一个点 ,得搞程见表 5-2,试拟合一曲面 ,确定合适的模型, 并由此找出最高点和该点的高程 .表 5-2 高程数据y x 100 200 300 40010020030040063669868066269771267462662463059855247847841233420x,ymeshgrid(100:100:400); z=636 698 680 662;697 712 674 626;624 630 598 552;478 478 412 334; mesh(x,y,z)hidden off
