-
Notifications
You must be signed in to change notification settings - Fork 188
/
go2.m
139 lines (125 loc) · 3.07 KB
/
go2.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
function [E,D] = go2(A,b)
%
% [E,D] = go2(A,b)
%
% Phase portrait of 2D linear o.d.e.
% u' = A * u + b
%
% b = 0 if not given
%
% E - eigenvectors of A
% D - eigenvalues
%
% Use cursor and mouse button to input initial data
% keyboard for other options:
%
% 'i' - type initial data
% 'e' or 'v' - plot eigendirections
% 'k' - keyboard
% 'p' - print eigenvalues and vectors
% 'q' - quit
% 'r' - reset zoom
% 't' - change time interval
% 'w' - zoom out
% 'z' - zoom in
% all others - print options
%
% See also LO2, GO3
colors = ['b','g','r','c','m','y'];
ncolors = size(colors,2);
ic = 1;
global Ao bo
if nargin == 1, bo = [0;0]; else bo = b; end
Ao = A;
[E,D] = eig(A);
if isreal(E)
e1 = D(1,1); e2 = D(2,2);
ev1 = E(:,1)'; ev2 = E(:,2)';
nv1 = ev1; nv2 = ev2;
realev = 1;
else
e1 = D(1,1);
ev = E(:,1)'; ev1 = real(ev); ev2 = imag(ev);
nv1 = ev1/norm(ev1); nv2 = ev2/norm(ev2);
realev = 0;
end
tso = 20; ts = [0 tso];
xfa = 10; xfo = 20; vfo = .04;
ax = xfa * [-1 1 -1 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);
idfl = 1;
else
disp('Initial data must be row vector with 2 entries');
idfl = 0;
end
else
idfl = 1;
end
if idfl == 1
u0 = [x0 y0]';
[t,y] = ode45('lo2',ts,u0);
hold on; plot(y(:,1),y(:,2),colors(ic)); hold off;
[t,y] = ode45('lo2',-ts,u0);
hold on; plot(y(:,1),y(:,2),colors(ic)); hold off;
xx = get(gax,'Xlim'); vf = vfo * xx(2);
v0 = Ao * u0 - bo; v0 = vf * v0/norm(v0);
wd = v0(1) - v0(2); ws = v0(1) + v0(2);
hold on;
plot([x0 - wd,x0,x0 - ws],[y0 - ws,y0,y0 + wd],colors(ic));
hold off;
ic = ic + 1; if ic > ncolors, ic = 1; end
end
case {101,118} % 'e' or 'v' - plot eigendirections
ez1 = xfo * nv1' *[-1 1]; ez2 = xfo * nv2' *[-1 1];
hold on; plot(ez1(1,:),ez1(2,:),'k');
plot(ez2(1,:),ez2(2,:),'k'); hold off;
case 107 % 'k' - keyboard
keyboard
case 112 % 'p' - print eigenvalues and vectors
if realev == 1
disp(['Eigenvalue: ', num2str(e1),' Eigenvector: [',num2str(ev1),']''']);
disp(['Eigenvalue: ', num2str(e2),' Eigenvector: [',num2str(ev2),']''']);
else
disp(['Eigenvalue: ', num2str(e1),...
' Eigenvector: [',num2str(ev1),']'' + i * [',num2str(ev2),']''']);
end
case 113 % 'q' - quit
stopit = 1;
case 114 % 'r' - reset zoom
zoom out; %vf = vfo;
case 116 % 't' - change time interval
tso = input('Input time interval :'); ts = [0 tso];
case 119 % 'w' - zoom out
zoom(.5); %vf = vf / .5;
case 122 % 'z' - zoom in
zoom(2.); %vf = vf / 2.;
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