下面是绘制图5.8的参考程序。 【 clf, N=linspace(0,300,50);
dN=0.3134*(1-N/250).*N; plot(N,dN),hold on, plot([0 300],[0,0]),
plot([0,250/2,250],[0,0,0],'o'),
xlabel('N','fontsize',11),ylabel('dN','fontsize',11),
text(N(32),dN(32),'\\leftarrow\\it{d N / d t}>0,相点递增右移','fontsize',11), text(125,dN(45),'\\it{dN/dt}<0,相点递减左移 \\rightarrow','fontsize',11);
h=text(251/2,1.5,'\\it{N_m/2}');set(h,'fontsize',11) 】
5.观察程序及其说明
zxy5_1.m (绘制函数图象，图5.3)
【 clf, x=sym('x'); f=(x-3)^2/(4*(x-1)); g=x/4-5/4;
hold on,
h=line([-8 8],[0,0]); set(h,'color’,'red’); h=line([0 0],[-8,8]); set(h,'color’,'red’); line([1 1],[-8 8]);plot([-1 1 3],[-2,0,0],'o’),
ezplot(g,[-8 8]); ezplot(f,[-8,8]), %符号函数绘图 text(-1-0.5,-2-0.5,'(-1,-2)’);text(1,0-0.5,'(1,0)'); text(3,0.5,'(3,0)'); x=1.4;text(x,subs(f),'\\leftarrow{(x-3)}^{2}/{4(x-1)}'); x=0.6;text(x,subs(f),'\\leftarrow{(x-3)}^{2}/{4(x-1)}'); x=2.5;text(x,subs(g),'\\leftarrow斜渐近线{y=x/4-5}');
text(1,-2,'\\leftarrow垂直渐近线x=1');title('(x-3)^2/4(x-1)') 】 zxy5_2.m （人口预测，图5.5） 【 global p1;clf,
t1=[1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 ];
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.0 72.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5]; p1(1)=3.9; p1(2)=250;p1(3)=1790;p1(4)=0.03134;
26
[t,N]=ode23('Logistic_fun',[1790 2100],3.9);
plot(t,N,t1,x,'o',t,250*ones(1,length(t))),axis([1790 2100 0 300]),
xlabel('t','fontsize',11), ylabel('N','fontsize',11) 】 Logistic_fun.m
【 function dN = Logistic_fun(t,N)
global p1
N0=p1(1);Nm=p1(2);t0=p1(3);r=p1(4); dN =r*(1-N/Nm)*N; 】
zxy5_3.m （方程的解轨线和相轨线，图5.6） 【 clear,clf
global p1;
p1(1)=3.9; p1(2)=250;p1(3)=0;p1(4)=0.03134;Nm=p1(2); tpas=linspace(0,300,1000);
plot([0 0],[0,500],':',[0 300],[Nm,Nm],':',[0 300],[Nm/2,Nm/2],':'), axis([-50 300 0 500]),xlabel('t'),ylabel('N'),hold on text(-30,Nm,'\\it{N_m}\\rightarrow'); text(-35,Nm/2,'\\it{N_m/2}\\rightarrow'); button=1; while button==1 k=[];
[t0,N0,button]=ginput(1);
[t,N]=ode23('Logistic_fun',tpas,N0); k=find(N<=Nm/2+1&N>=Nm/2-1);
ts=tpas(k);Ns=N(k); text(-35,N0,'\\it{N_0}\\rightarrow'); plot(t0,N0,'o',ts,Ns,'square',t,N,':'),hold on comet(t,N),pause,comet(0*ones(length(t),1),N) end 】 这一程序是不难读懂的。
5.3 应用、思考和练习
27
5.3.1.函数作图
◆（2） 下面的绘图较复杂一些，是一个很好的练习。 zxy5_4.m
【 clf,n=2000;a=-4;b=6;c=-8;d=8;
t=linspace(a,b,n);
x=(t.^2)./(t-1);y=t./(t.^2-1); kx=find(abs(x)>=d);x(kx)=NaN; ky=find(abs(y)>=d);x(ky)=NaN; plot(t,[x;y],'.','markersize',3), hold
on,plot([a
b],[0,0],'r',[0
0],[c,d],':'),axis([a b c d]), xlabel('t'),ylabel('x and y') ,
text(-3.8,7,'put any key to show x=x(t)');pause,comet(t,x),
text(-3.8,6,'put any key to show y=y(t)');pause,comet(t,y) 】
图5.9 x(t),y(t)的曲线
5.3.2.平衡点的分类 5.3.3定性分析的应用
1.捕鱼业持续的收获
画定性分析图的程序zxy5_4.m
【 clf,clear,N=50; r=4.4;E=[0.5 2.2 6.5];x=linspace(0,N,30);
f1=r*x.*(1-x/N);plot(x,r.*x,':','linewidth',2),axis([0 50 0 80]),hold on text(x(10),r*x(10),['\\leftarrow y=rx, r= ',num2str(r)]) for i=1:3 f2(i,:)=E(i)*x;
text(x(5),f2(i,5),['\\leftarrow y=',num2str(E(i)),'x']) end
28
plot(x,f1,x,f2),hold on, text(x(22),f1(22),['\\leftarrow
dx/dt=rx(1-x/N)']),xlabel('x'),ylabel('dx/dt') 】
2. 蚜虫生长和跃变
实验６．最优化实验
6.1知识要点与背景
6.1.1 由简入繁： 最佳水槽断面问题的推广
6.1.2 微分法求最大和最小
◆问题4 （1）求受检点： zxy6_1.m
【 syms x1 x2 %
定义符号变量。
?f2?3x1?6x1?9， ?x1f=x1^3-x2^3+3*x1^2+3*x2^2-9*x1; % 函数z。
v=[x1 x2];df=jacobian(f,v) %计算雅可比
?(z)/?(x1,x2)??z。
[X,Y]=solve(df(1),df(2));[X,Y] % 用指令solve求驻
点。 】
zxy6_2.m
【 clf,xmin=-5;xmax=3.5;ymin=-3;ymax=5;
x1=linspace(xmin,xmax,30);x2=linspace(ymin,ymax,30);
[X1,X2]=meshgrid(x1,x2);Z= X1.^3.-X2.^3+3*X1.^2+3*X2.^2-9*X1; contour(X1,X2,Z,60);,hold on, xp=[-3,1,-3,1];yp=[0 0 2 2]; plot(xp,yp,'ro'),axis([xmin xmax ymin ymax]),colorbar
29
xlabel('x_1'),ylabel('x_2'), for i=1:length(xp)
text(xp(i),yp(i),['\\leftarrow (',num2str(xp(i)),',',num2str(yp(i)),')'] ) end 】
6.2 实验与观察(Ⅰ)：模拟盲人下山的迭代寻优法
zxy6_3.m（盲人下山的模拟） 【 clf, 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]=zxy6_3f(X1,X2); %计算函数和梯度向量。 contour(X1,X2,Z,30), %画等值线图。 axis([xmin xmax ymin ymax]),hold on,
axis equal, %该命令将使横轴、纵轴
具有相同比例，避免失真。
plot([1.46808510638298],[1.148936170212776],'o'), %标注最优点。 axis([xmin xmax ymin ymax])
x=[];y=[]; %开始用鼠标选点，按左键选点，按右键中止选点过
程。
disp('Select a point by put on mouse left-key')
%disp指令,在命令窗口显示文字。 disp('Stop selecting point by put on mouse right-key')
button=1; %button和ginput命令结合使用可用鼠标选点， 按左键时
button＝1。
x=[];y=[];
while button==1 [xi,yi,button]=ginput(1);
%ginput(n)用鼠标选n个点，xi,yi分别为点的横坐标
和纵坐标。
plot([xi],[yi],'r.','MarkerSize',10)，hold on, %画所选的
30