1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
function []=portrait(fch,odem,xdim,npts,farrow,pinit)
// phase portrait
// Copyright INRIA
xselect();
ncnl=lines();
lines(0);
//addtitle(fch);
[lhs,rhs]=argn(0);
// minimal calling sequence
if rhs<=1,odem='default';end
// Interactive version
if odem == 'discrete'; style_d=x_choose(['trait continu','points'],['option de dessin']);
style_d=mini(-style_d,1);end
if rhs <= 2,
if ~isdef('p_xdim');p_xdim=['-1';'-1';'1';'1'];end
rep=x_mdialog('Graphic boundaries',...
['xmin';'ymin';'xmax';'ymax'],...
p_xdim);
if rep<>[] ;p_xdim=rep;end
xdim=evstr(p_xdim');
// Test sur le cadre
if xdim(3) <= xdim(1),
write(%io(2),'Error: xmin < xmax '),lines(ncnl(1));return,end
if xdim(4) <= xdim(2),
write(%io(2),'Error: ymin < ymax '),lines(ncnl(1));return,end
end
res=x_choose(['yes';'no'],'Do you also want to draw the vector field')
if res==1;
if ~isdef('p_nxx');p_nxx=['10';'10'];end;
rep=x_mdialog('Number of grid points',...
['Nx';'Ny'],p_nxx);
if rep<>[] then p_nxx=rep ;end
nxx=evstr(p_nxx);
nx=maxi(nxx(1),2)
ny=maxi(nxx(2),2)
nx=(0:(nx-1))/(nx-1)
ny=(0:(ny-1))/(ny-1)
nx=xdim(1)*(ones(nx)-nx)+xdim(3)*nx;
ny=xdim(2)*(ones(ny)-ny)+xdim(4)*ny;
fchamp(fch,0,nx,ny,1.0,xdim);
style="000";
else
p_nxx=['10';'10'];
style="051";
end
plot2d(0,0,0,style," ",xdim);
if rhs<=3,
if ~isdef('p_npts');p_npts=['100';'0.1'];end;
rep=x_mdialog('Requested points and step ',...
['n points';'step'],p_npts);
if rep <> [] then p_npts=rep;end
npts=evstr(p_npts');
end
ylast=(1/2)*[xdim(3)+xdim(1),xdim(4)+xdim(2)]';
if rhs<=4,farrow='f';end;
if rhs<=5
// Loop on Initial points
go_on=1
while go_on==1,
ftest=1;
while ftest==1,
n=x_choose(['New initial point';'Continue ode';'Quit'],"Choose ");
n=n-1;
if n==-1,go_on=0;lines(ncnl(1));
[p_xdim,p_npts,p_nxx]=resume(p_xdim,p_npts,p_nxx);return;end
if n==2,go_on=0;lines(ncnl(1));
[p_xdim,p_npts,p_nxx]=resume(p_xdim,p_npts,p_nxx);return;end
if n==0,[i,x,y]=xclick(); x0=[x,y];end;
if n==1,x0=ylast';end;
ftest=desorb(odem,x0',npts,fch,farrow,xdim);
if ftest==1;x_message('Initial value out of boundaries'),end
end
end
else
// No question mode
res=desorb(odem,pinit,npts,fch,farrow,xdim);
if res==1,write(%io(2),'Points hors du cadre elimines ');end;
end
lines(ncnl(1));
[p_xdim,p_npts,p_nxx]=resume(p_xdim,p_npts,p_nxx);
function []=addtitle(fch)
// Adds know titles
//!
if type(fch)<>11& type(fch)<>13 then return;end;
if fch==linear,xtitle("Systeme lineaire"," "," ",0);end
if fch==linper,xtitle("Systeme lineaire perturbe "," "," ",0);end
if fch==cycllim,xtitle("Systeme avec cycle limite "," "," ",0);end
if fch==bioreact,xtitle("Bioreacteur ","biomasse ","sucre ",0);end
if fch==lincom,xtitle("Systeme lineaire commande "," "," ",0);end
if fch==p_p,xtitle("Modele proie-predateur ","proies ","predateurs ",0);end
if fch==compet,xtitle("Modele de competition ","population 1 "...
,"population2 ",0);end
if fch=='bcomp',xtitle("Modele de competition observe-comtrole ",...
"population 1 ","population2 ",0);end
if fch=='lcomp',xtitle("Modele de competition linearise observe-comtrole ",...
"population 1 ","population2 ",0);end
function [res]=desorb(odem,x0,n1,fch,farrow,xdim);
// Used by portrait
//!
res=0
[nn1,n2]=size(x0);
style=1;
if odem=='discrete', style=style_d;end
for i=1:n2,
ftest=1;
if x0(1,i) > xdim(3), ftest=0;end
if x0(1,i) < xdim(1), ftest=0;end
if x0(2,i) > xdim(4), ftest=0;end
if x0(2,i) < xdim(2), ftest=0;end
if ftest==0;res=1,ylast=x0,else
write(%io(2),'Calling ode')
if odem=='default' then
y=ode([x0(1,i);x0(2,i)],0,n1(2)*(0:n1(1)),fch);
else
y=ode(odem,[x0(1,i);x0(2,i)],0,n1(2)*(0:n1(1)),fch);
end;
[nn1,n11]=size(y);
// on coupe la trajectoire au temps d'arret T
// T d'atteinte du bord du cadre
[mi1,ki1]=mini(y(1,:),xdim(3)*ones(1,n11));
[ma1,ka1]=maxi(y(1,:),xdim(1)*ones(1,n11));
k1=maxi(ki1,ka1);
[mi2,ki2]=mini(y(2,:),xdim(4)*ones(1,n11));
[ma2,ka2]=maxi(y(2,:),xdim(2)*ones(1,n11));
k2=maxi(ki2,ka2);
[m11,k11]=maxi(k1);
[m22,k22]=maxi(k2);
if k11==1,k11=n1(1);end
if k22==1,k22=n1(1);end
kf=mini(k11,k22);
if kf==1, kf=n1(1),end
if farrow=='t',
plot2d4("gnn",y(1,1:kf)',y(2,1:kf)',style,"000"," ",xdim);
else
plot2d(y(1,1:kf)',y(2,1:kf)',style,"000"," ",xdim);
end,
ylast=y(:,kf);
end
end
[ylast]=resume(ylast)
|
