【实验目的】

1．了解二元函数图形的制作。

2．空间曲面等高线的制作。

3．学习掌握MATLAB软件有关的命令。

【实验内容】

画出函数的图形,并画出其等高线。

【实验准备】

1．曲线绘图的MATLAB命令

MATLAB中主要用mesh,surf命令绘制二元函数图形。

mesh(x,y,z)  画网格曲面，这里x,y,z是三个数据矩阵，分别表示数据点的横坐标，纵坐标和函数值，该命令将数据点在空间中描出,并连成网格。

【实验方法与步骤】

练习1  画出函数的图形,不妨将区域在。用MATLAB作图的程序代码为：

>>clear;

>>x=-3:0.1:3;  %x的范围为[-3,3]

>>y=-3:0.1:3; %y的范围为[-3,3]

>>[X,Y]=meshgrid(x,y); %将向量x,y指定的区域转化为矩阵X,Y

>>Z=sqrt(X.^2+Y.^2);    %产生函数值Z

>>mesh(X,Y,Z)

结果如图5.1。图5.1是网格线图,如果要画完整的曲面图,只需将上述的MATLAB代码mesh(X,Y,Z)改为surf(X,Y,Z), 结果如图5.2

图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);

>>contour3(X,Y,Z,10)  %画10条等高线

>>xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis')  %三个坐标轴的标记

>>title('Contour3 of Surface')  %标题

>>grid on  %画网格线

结果如图5.3.

图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.

图5.4  等高线

如果要画的等高线，则用命令

>>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,[1 1])

结果如图5.5。

图5.5  等高线

练习1中，函数值可简单算出。在有些情况下，函数值不能简单算出。这是因为x和y的值可能是非均匀间隔的甚至是随机分布的，也可能使用了不同的坐标系，比如非长方形的网。出现这些情况时,MATLAB中的函数griddata就用来产生经查值后的均匀间隔数据以作图。

练习２ 二次曲面的方程如下

讨论参数对其形状的影响。

    本练习的关键在于如何作出三维曲面图形，特别注意在给定值求时，若有开方运算，一是会出现虚数，二是对实数也有正负两个解。为了使虚数不出现在绘图中，采用了一种技巧，就是将虚数都换成非数(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网格坐标

>>[x,y]=meshgrid(xgrid,ygrid); %确定个点的x,y网格坐标

>>z=c*sqrt(d-y.*y/b^2-x.*x/a^2); u=1; %u=1，表示z要取正值

>>z1=real(z);  %取z的实部z1

>>for k=2:N-1  %以下7行程序的作用是取消z中含虚数的点

     >>for j=2:N-1

        if imag(z(k,j))~=0 z1(k,j)=0; end

        if all(imag(z([k-1:k+1],[j-1:j+1])))~=0 za(k,j)=NaN; end

        end

     end

     >>surf(x,y,z1), hold on  %画空间曲面

     >>if u==1 z2=-z1; surf(x,y,z2);  %u=1时加画负半面

        axis([-abs(a),abs(a), -abs(b), abs(b), -abs(c), abs(c)]);

     end

     >>xlabel('x'),ylabel('y'),zlabel('z')

     >>hold off

运行程序，当时的结果见图5.6, 

当时的结果见图5.7,

当时的结果见图5.8,

图5.6  椭球面

图5.7  双曲面

图5.8  椭球双曲面

练习3 列出求空间两任意曲面的交线的程序。

两空间曲面方程连立起来，就形成一个空间曲线的方程。这个曲线能满足两个曲面的方程，因而也就是这两个空间曲面的交线。显示这两个曲面并不难，用两次mesh语句即可,但要显示其交线，必须先找到各个交点，因为数值计算得到的是离散点，难以找到两个曲面上完全重合的点，本程序采用了设置限的方法，只要在同一网格点处，两曲面的z之之差小于设定限，就认为它是交点，限值设定几次要才能定的好。

     下面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); hold; mesh(x,y,z2);  %再一个图上同时画出两个曲面

>>r0=(abs(z1-z2)<=0.1);  %求两曲面z坐标差小于0.1的网格矩阵

>>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  两曲面的交线

如果想改表曲面方程，可以在程序中改动第二行和第三行。但这样的程序还不是通用的，最好程序运行时能向用户提问，允许用户输入曲面方程。此时就要用到字符串功能和eval命令。

s1=input(‘输入第一个方程’,’s’);

在原来的z1方程语句处改为z1=eval(s1);类似地输入第二个方程。此外，应使用户能给出定义域和间隔。这实现起来比较简单，只要把第一句改为

[x,y]=meshgrid(xmin:dx:xmax,ymin:dy:ymax);

其中，xmin，dx，xmax，ymin，dy，ymax可由程序给出屏幕提问，让用户用键盘输入。当然，这样又增加了运行时的麻烦，所以编程时要找一个折衷的选择，要有一定的灵活性又不能太麻烦，应恰到好处。

练习4 用平行界面法讨论由方程构成的马鞍面形状。

我们只需对练习3种的程序作如下修改：

定义域网格改为[x,y]=meshgrid(-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;  %第一个曲面方程

>>a=input('a=(-50>>subplot(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)<=1;  %求两曲面z坐标差小于1的网格

>>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 %使得第二个分图取第一个分图的坐标系

    执行此程序,并输入a=8,得到的三维图形及交线见图5.10, 当a=-20,得到的三维图形及交线见图5.11,可见从上而下,其横切面交线发生了很大的变化.

图5.10  马鞍面的水平截面(a=8)

图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程序代码为：

>>clear;

>>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.

图5.12  柄图

显然上面数据是残缺不全的，下面用插值的方法画出完整的曲面，相应的MATLAB程序代码为：

>>xi=0:0.2:3; yi=2:0.2:4;   %选定x,y的范围

>>[X,Y]=meshgrid(xi,yi); %产生网格向量X,Y

>>Z=griddata(x,y,z,X,Y,'cubic');    %’cubic’采用三角形三次插值

>>mesh(X,Y,Z); title('Griddata');

>>xlabel('x'),ylabel('y'),zlabel('z')

结果如图5.13.

图5.13  插值曲面

练习6 （海底测量）表5-1给出水面直角坐标(x,y)处水深z,这时在低潮时测得的。如果船的吃水深度为5米，试问在矩形域中船应避免进入那些区域？

表5-1  水深数据

x(m)

y(m)

>>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.

图5.14 测量点的位置

由图5.8可见，这是一批不规则数据。由于没有先验函数，我们使用插值法。为了使结果更直观，考虑将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米的两个危险区域：

>>contour(X,Y,H, [-5,-5],'k')  %’k’表示等高线的颜色为黑色

图5.16 两个危险区域 

【练习与思考】

1.画出空间曲线在范围内的图形，并画出相应的等高线。

>>[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）椭球面

>> [x,y,z]=ellipsoid(0,0,0,3,2,1,100);

>> mesh(x,y,z)

>>axis equal

（2）椭圆抛物面

>> syms u v

>> ezmesh(3*u*sin(v),2*u*cos(v),4*u^2)

（3）单叶双曲面

>> [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)

（4）双曲抛物面

>> [x,y]=meshgrid(-10:0.1:10);

>> z=(x.^2-y.^2)/3;

>> mesh(x,y,z)

（5）旋转面

>>[t,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)

（6）圆锥面

>> [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）环面

>> syms u v

>> x=(3+0.4*cos(u))*sin(v);

>> y=(3+0.4*cos(u))*cos(v);

>> z=sin(u)*3;

>> z=0.4*sin(u);

>> ezmesh(x,y,z)

（8）正螺面

>> 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  高程数据

>> z=[636 698 680 662;697 712 674 626;624 630 598 552;478 478 412 334];

>> mesh(x,y,z)

>>hidden off下载本文