matlab画图.doc

实验五 二元函数的图形
【实验目的】
1.了解二元函数图形的绘制。
2.了解空间曲面等高线的绘制。
3.了解多元函数插值的方法。
4.学习、掌握MATLAB软件有关的命令。
【实验内容】
画出函数的图形，并画出其等高线。
【实验准备】
1.曲线绘图的MATLAB命令
MATLAB中主要用mesh,surf命令绘制二元函数图形。主要命令
mesh（x，y，z）画网格曲面，这里x，y，z是数据矩阵，分别表示数据点的横坐标，纵坐标和函数值，该命令将数据点在空间中描出，并连成网格。
surf（x，y，z）画完整曲面，这里x，y，z是数据矩阵，分别表示数据点的横坐标，纵坐标和函数值，该命令将数据点所表示曲面画出。
【实验重点】
1. 二元函数图形的描点法
2. 曲面交线的计算
3. 地形图的生成
【实验难点】
1. 二元函数图形的描点法
2. 曲面交线的计算
【实验方法与步骤】
练习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.3
如果画等高线，用contour，contour3命令。
contour画二维等高线。
contour3画三维等高线。画图5.3所示的三维等高线的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.4

如果画图5.4所示的二维等高线，相应的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('Contour3 of Surface')
>>grid on
运行结果为

如果要画z=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,[1 1])
运行结果为

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

讨论参数a,b,c对其形状的影响。
相应的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); %确定N×N个点的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 z1(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
运行程序，当a=5,b=4,c=3,d=1,N=50时结果为

当a=5i,b=4,c=3,d=1,N=15时结果为

当a=5i,b=4i,c=3,d=0.1,N=10时结果为

【练习与思考】
1. 画出空间曲面在范围内的图形，并画出相应的等高线。
解：
clear;close;
u=-30:0.5:30;
v=-30:0.5:30;
[x,y]=meshgrid(u,v);
z=10*sin(sqrt(x.^2+y.^2))./sqrt(1+x.^2+y.^2);
subplot(1,2,1);
mesh(x,y,z)
subplot(1,2,2);
contour(x,y,z,10)

2. 根据给定的参数方程，绘制下列曲面的图形。
a) 椭球面，，；
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=3.*cos(U).*sin(V);
Y=2.*cos(U).*cos(V);
Z=sin(U);
surf(X,Y,Z);
axis equal

b) 椭圆抛物面，，；
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=3.*U.*sin(V);
Y=2.*U.*cos(V);
Z=4.*U.^2;
surf(X,Y,Z);
shading interp;
colormap(hot);
axis equal

c) 单叶双曲面，，；
解：
clear;close;
u=-2:0.1:2;
v=-2:0.1:2;
[U,V]=meshgrid(u,v);
X=3*sec(U).*sin(V);
Y=2*sec(U).*cos(V);
Z=4*tan(U);
mesh(X,Y,Z);
shading interp;
colormap(jet);
axis equal

d) 双叶抛物面，，；
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=U;
Y=V;
Z=(U.^2-V.^2)/3;
mesh(X,Y,Z);
shading interp;
colormap(jet);
axis equal

e) 旋转面，，；
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=log(U).*sin(V);
Y=log(U).*cos(V);
Z=U;
mesh(X,Y,Z);
shading interp;
colormap(jet);
axis equal

f) 圆锥面，，；
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=U.*sin(V);
Y=U.*cos(V);
Z=U;
mesh(X,Y,Z);
shading interp;
colormap(jet);
axis equal

g) 环面，，；
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=(3+0.4*cos(U)).*cos(V);
Y=(3+0.4*cos(U)).*sin(V);
Z=0.4*sin(V);
mesh(X,Y,Z);
shading interp;
colormap(jet);
axis equal

h) 正螺面，，。
解：
clear;close;
u=-4:0.1:4;
v=-4:0.1:4;
[U,V]=meshgrid(u,v);
X=U.*sin(V);
Y=U.*cos(V);
Z=4*V;
mesh(X,Y,Z);
shading interp;
colormap(jet);
axis equal

3. 在一丘陵地带测量高程，和方向每隔100米测一个点，得高程见表5-2，试拟合一曲面，确定合适的模型，并由此找出最高点和该点的高程。
表5-2 高程数据
y x
100
200
300
400
100
200
300
400
636
698
680
662
697
712
674
626
624
630
598
552
478
478
412
334
解：
clear;close;
x=[100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400];
y=[100 200 300 400 100 200 300 400 100 200 300 400 100 200 300 400];
z=[636 697 624 478 698 712 630 478 680 674 598 412 662 626 552 334];
xi=100:5:400;
yi=100:5:400;
[X,Y]=meshgrid(xi,yi);
H=griddata(x,y,z,X,Y,'cubic');
surf(X,Y,H);
view(-112,26);
hold on;
maxh=vpa(max(max(H)),6)
[r,c]=find(H>=single(maxh));
stem3(X(r,c),Y(r,c),maxh,'fill')