实验五 二元函数的图形
【实验目的】
1. 了解二元函数图形的制作。
2. 空间曲面等高线的制作。
3. 了解多元函数插值的方法。
4. 学习掌握MATLAB软件有关的命令。
【实验内容】
画出函数
的图形,并画出其等高线。
【实验准备】
1.曲线绘图的MATLAB命令
MATLAB中主要用mesh,surf命令绘制二元函数图形。
mesh(x,y,z) 画网格曲面,这里x,y,z是三个数据矩阵,分别表示数据点的横坐标,纵坐标和函数值,该命令将数据点在空间中描出,并连成网格。
surf(x,y,z) 画完整曲面,这里x,y,z是三个数据矩阵,分别表示数据点的横坐标,纵坐标和函数值,该命令将数据点所表示曲面画出。
可以用help mesh, help surf查阅有关这些命令的详细信息
【实验方法与步骤】
练习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就用来产生经查值后的均匀间隔数据以作图。
练习2 二次曲面的方程如下
讨论参数
对其形状的影响。
本练习的关键在于如何作出三维曲面图形,特别注意在给定
值求
时,若有开方运算,一是会出现虚数,二是对实数也有正负两个解。为了使虚数不出现在绘图中,采用了一种技巧,就是将虚数都换成非数(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)
z(m)
129
7
4
140
141
8
108
28
6
88
147
8
185
22
6
195
137
8
105
85
8
x(m)
y(m)
z(m)
157
-6
9
107
-81
9
77
3
8
145
45
8
162
-66
9
162
84
4
117
-38
9
我们首先看测量点的位置:
>>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. 画出空间曲线
在
范围内的图形,并画出相应的等高线。
2. 根据给定的参数方程,绘制下列曲面的图形。
(1) 椭球面
(2) 椭圆抛物面
(3) 单叶双曲面
(4) 双曲抛物面
(5) 旋转面
(6) 圆锥面
(7) 环面
(8) 正螺面
3. 在一丘陵地带测量搞程,x和y方向每隔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
PAGE
2
_1124691745.unknown
_1131476341.unknown
_1131477141.unknown
_1131477377.unknown
_1131477890.unknown
_1131511815.unknown
_1131477376.unknown
_1131476544.unknown
_1131476553.unknown
_1131476409.unknown
_1124695309.unknown
_1124695369.unknown
_1124695475.unknown
_1124695509.unknown
_1124695596.unknown
_1124695423.unknown
_1124695339.unknown
_1124694492.unknown
_1124695266.unknown
_1124694450.unknown
_1124628874.unknown
_1124688433.unknown
_1124688741.unknown
_1124628888.unknown
_1124628855.unknown