-
Notifications
You must be signed in to change notification settings - Fork 192
/
phaseg.m
123 lines (105 loc) · 2.29 KB
/
phaseg.m
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
function phaseg(o,ti)
%
% phasem(o,ti)
%
% Phase portrait of 2D o.d.e. given in file o
% ti - time interval
%
% Use cursor and mouse button to input initial data
% keyboard for other options:
%
% 'i' - type initial data
% 'k' - keyboard
% 'q' - quit
% 'r' - reset zoom
% 't' - change time interval
% 'w' - zoom out
% 'z' - zoom in
% all others - print options
%
% See also GO2
if nargin < 2, ti = 50; end
pa = 10;
maxy = 1;
colors = ['b','g','r','c','m','y'];
ncolors = size(colors,2);
ic = 1;
tso = ti; ts = [0 tso];
axo = 4 * [-1 1 -1 1];
zfac = 1;
ax = zfac * axo;
vfo = .1;
figure(2); clf;
figure(1);clf;
gax = axes;
axis(ax);
zoom on;
stopit = 0;
while stopit == 0
[x0 y0 n] = ginput(1);
if size(n) == [0 0]
stopit = 1;
else
switch n
case {1,105} % mouse button - input initial data
% 'i' - keyboard initial data
if n == 105
id = input('Initial Data : ');
if size(id) == [1 2]
x0 = id(1); y0 = id(2);
disp(['Initial data : ',num2str(x0),' ',num2str(y0)]);
idfl = 1;
else
disp('Initial data must be row vector with 2 entries');
idfl = 0;
end
else
disp(['Initial data : ',num2str(x0),' ',num2str(y0)]);
idfl = 1;
end
if idfl == 1
u0 = [x0 y0]';
[t,y] = ode45(o,ts,u0);
hold on; plot(y(:,1),y(:,2),colors(ic)); hold off;
% [t,y] = ode45(o,-ts,u0);
% hold on; plot(y(:,1),y(:,2),colors(ic)); hold off;
figure(2);
maxi = max(abs(y(:,1))) + .5;
maxy = max(maxy,maxi);
hold on; plot(t,y(:,1),colors(ic)); hold off;
axis([0,tso,-maxy,maxy]);
figure(1);
ic = ic + 1; if ic > ncolors, ic = 1; end
end
case 107 % 'k' - keyboard
keyboard
case 113 % 'q' - quit
stopit = 1;
case 114 % 'r' - reset zoom
zfac = 1;
ax = zfac * axo;
axis(ax);
case 116 % 't' - change time interval
tso = input('Input time interval :'); ts = [0 tso];
case 119 % 'w' - zoom out
zfac = 2* zfac;
ax = [x0 x0 y0 y0 ] + zfac * axo;
axis(ax);
case 122 % 'z' - zoom in
zfac = zfac/2;
ax = [x0 x0 y0 y0 ] + zfac * axo;
axis(ax);
otherwise
disp( 'Options are:')
disp('i - type initial data')
disp('e or v - plot eigendirections')
disp('k - keyboard')
disp('t - change time interval')
disp('p - print eigenvalues and vectors')
disp('r - reset zoom')
disp('w - zoom out')
disp('z - zoom in')
disp('q - quit')
end
end
end