1701009615
编写MATLAB程序绘制函数如下：
1701009616
1701009617
clc,clear,close all %清屏和清除变量 warning off %消除警告 %盲人下山 figure(‘color’,[1,1,1]) %设置图形背景为白色 ezmesh(‘8*x1^2+9*x2^2-8*x1*x2-12*x1-6*x2’,[-2,4,-2,4]) %曲面绘制 view([-24.5 2]); %图形视角 grid off %不显示网格 axis off %不显示坐标轴
1701009618
1701009619
运行程序输出图形如图8-3所示。
1701009620
1701009621
1701009622
1701009623
1701009624
图8-3 模拟山峰
1701009625
1701009626
如图8-3所示，对于该函数图像，可以等效为很多个等值线。盲人下山问题即可简化为盲人在不同的等值线上的运动轨迹，为了较为清晰地模拟盲人下山的随机性，采用鼠标随机选点，大家可以任意选取属于自己的下山路径，具体的编程如下：
1701009627
1701009628
clc,clear,close all %清屏和清除变量 warning off %消除警告 %盲人下山 a=-2;b=4; xmin=a;xmax=b;ymin=a;ymax=b; x1=linspace(xmin,xmax,100); %等分区间 x2=linspace(ymin,ymax,100); %等分区间 [X1,X2]=meshgrid(x1,x2); %数据栅格化 [Z,DZ1,DZ2]=jt(X1,X2); %调用函数 axis([xmin xmax ymin ymax]); axis equal plot([1.46808510638298],[1.148936170212776],‘o’) %画图 axis([xmin xmax ymin ymax]); %坐标轴设置 hold on %图形保持句柄 disp(‘鼠标左键选取等高线’) %显示字符串 disp(‘停止下山’) %显示字符串 button=1; x=[];y=[]; while button==1 [xi,yi,button]=ginput(1); %鼠标选点 plot([xi],[yi],‘r.’,‘markersize’,10); %画点图 hold on; %图形保持句柄 [zi,dz1,dz2]=jt(xi,yi); %调用函数 v=zi; contour(X1,X2,Z,[v v],’-‘); %等高线图 axis([xmin xmax ymin ymax]); %坐标轴设置 x=[x,xi];y=[y,yi]; H_line2=plot(x,y); set(H_line2,‘color’,‘red’,‘linewidth’,2); %设置线宽和线的颜色 xt=xi-dz1;yt=yi-dz2; H_line=plot([xi xt],[yi yt],‘k:’,‘linewidth’,1); %画图 end
1701009629
1701009630
运行程序输出图形如图8-4所示。
1701009631
1701009632
1701009633
1701009634
1701009635
图8-4 盲人下山的路径
1701009636
1701009637
盲人起点如图8-4所示中的小圆圈点，图8-4中的线表示等高度线，红色的粗实线表示盲人选择的下山路径，当然也可以选择如图8-5和图8-6所示的下山路径。
1701009638
1701009639
1701009640
当然对于这个函数而言，给定一个初始点（x，y）和一个目标点（0，0），是存在一条最优路径的。
1701009641
1701009642
1701009643
1701009644
图8-5 下山路径选择1 图8-6 下山路径选择2 绘制该函数的梯度图像，程序如下：
1701009645
1701009646
xmin=0; %x最小值设置 xmax=4; %x最大值设置 ymin=0; %y最小值设置 ymax=4; %y最大值设置 x1=linspace(xmin,xmax,20); %等分区间 x2=linspace(ymin,ymax,20); %等分区间 [X1,X2]=meshgrid(x1,x2); %数据栅格化 [Z,DZ1,DZ2]=jt(X1,X2); %调用函数 figure(‘color’,[1,1,1]) %设置背景为白色 contour(X1,X2,Z,40) %等高线图 hold on %图形保持句柄 quiver(X1,X2,-DZ1,-DZ2) %绘制箭头
1701009647
1701009648
运行程序输出图形如图8-7所示。
1701009649
1701009650
1701009651
1701009652
1701009653
图8-7 等高线和曲面的负梯度图
1701009654
1701009655
从图8-7中可清晰的知道，每一个起点位置到山顶峰点之间的最快梯度路径，同样，对于该函数很容易求出一条最优路径，选取起点为（0.9,1.2,-9.63），终点为（0,0,0），编程如下：
1701009656
1701009657
clc,clear,close all %清屏和清除变量 warning off %消除警告 %盲人下山 [X1,Y1]= meshgrid(-2:0.1:4,-2:0.1:4); %数据平面栅格化 x1 = X1; x2 = Y1; Z1 = 8*x1.^2+9*x2.^2-8.*x1.*x2-12*x1-6*x2; %函数 figure(‘color’,[1,1,1]) %设置图形背景为白色 mesh(X1,Y1,Z1) %3D曲面绘制 %%NxM 矩阵化为一行 nX1 = size(X1); kk=1; for i=1:nX1(1,1) for j=1:nX1(1,2) X2(1,kk) = X1(i,j); %将矩阵化成一列 Y2(1,kk) = Y1(i,j); Z2(1,kk) = Z1(i,j); kk=kk+1; end end %%距离 for i=1:length(X2) for j=1:length(X2) Distance(i,j)=sqrt( (X2(i)-X2(j))^2 + (Y2(i)-Y2(j))^2 + (Z2(i)- Z2(j))^2 ); Distance(j,i)=Distance(i,j); end end %%邻接矩阵 A=zeros(length(X2),length(X2)); %构建邻接矩阵 %第一行数据 for i = 1:nX1(1,2)-1 A(i,i+1)=1; A(i,nX1(1,2)+i)=1; A(i,nX1(1,2)+i+1)=1; end for i=2:nX1(1,1)-1 %第二行到倒数第二行 for j=2:nX1(1,2)-1 %倒数第二列 A( (i-1)*nX1(1,2)+j, (i-2)*nX1(1,2)+j ) = 1; %该点正上方点 A( (i-1)*nX1(1,2)+j, (i-2)*nX1(1,2)+j-1 ) = 1; %该点正上方点，向左移一位 A( (i-1)*nX1(1,2)+j, (i-2)*nX1(1,2)+j+1 ) = 1; %该点正上方点，向右移一位 A( (i-1)*nX1(1,2)+j, (i-1)*nX1(1,2)+j-1 ) = 1; %该点行上，向左移一位 A( (i-1)*nX1(1,2)+j, (i-1)*nX1(1,2)+j+1 ) = 1; %该点行上，向右移一位 A( (i-1)*nX1(1,2)+j, (i)*nX1(1,2)+j ) = 1; %该点正下方点 A( (i-1)*nX1(1,2)+j, (i)*nX1(1,2)+j-1 ) = 1; %该点正下方点，向左移一位 A( (i-1)*nX1(1,2)+j, (i)*nX1(1,2)+j+1 ) = 1; %该点正下方点，向右移一位 end end for i=1:length(X2) for j=1:length(X2) A(j,i)=A(i,j); end end %%求最短路径 %起点（0.9,1.2,-9.63），终点（0,0,0） %画图 figure(‘color’,[1,1,1]) mesh(X1,Y1,Z1) hold on plot3(0.9,1.2,-9.63,‘r.’,‘Markersize’,40) plot3(0,0,0,‘g.’,‘Markersize’,40) hold off %% %起点在30行33列，终点在21行21列 D=Distance.*A; %相连节点计算距离 D(find(D==0))=99999; %两点无边相连时赋值为inf [Q_path, dmin] = dijkstra( 29*nX1(1,2)+33,20*nX1(1,2)+21 ,D); disp([‘最短路径为：’, num2str(dmin)]) %% figure(‘color’,[1,1,1]) mesh(X1,Y1,Z1) xlabel(‘X’);ylabel(‘Y’);zlabel(‘Z’) hold on plot3(0.9,1.2,-9.63,‘b.’,‘Markersize’,40) plot3(0,0,0,‘b.’,‘Markersize’,40) for i=1:length(Q_path) row = ( Q_path(i) - mod(Q_path(i), nX1(1,2)) )/nX1(1,2); column = mod(Q_path(i),nX1(1,2)); X3(i) = X1(column,row); Y3(i) = Y1(column,row); Z3(i) = Z1(column,row); plot3(X3(i),Y3(i),Z3(i),‘r.’,‘Markersize’,40,‘linewidth’,3); end %% figure(‘color’,[1,1,1]) mesh(X1,Y1,Z1) xlabel(‘X’);ylabel(‘Y’);zlabel(‘Z’) hold on plot3(0.9,1.2,-9.63,‘r.’,‘Markersize’,40) plot3(0,0,0,‘g.’,‘Markersize’,40) plot3(X3,Y3,Z3,‘b-‘,‘Markersize’,40,‘linewidth’,3); set(gca,‘Xdir’,‘reverse’); set(gca,‘Ydir’,‘reverse’);
1701009658
1701009659
运行程序输出结果如下：
1701009660
1701009661
最短路径为：9.8305
1701009662
1701009663
输出图形如图8-8所示。
1701009664