-
Notifications
You must be signed in to change notification settings - Fork 192
/
go3.m
79 lines (65 loc) · 1.32 KB
/
go3.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
function [E,D] = go3(A,dt,et)
%
% [E,D] = go3(A,dt,et)
%
% Phase portrait of 3D linear o.d.e.
% u' = A * u
%
% Optional arguments:
% dt - forward time interval
% et - backward time interval
%
% E - eigenvectors of A
% D - eigenvalues
%
% initial data input is either a
% k by 3 matrix giving k initial points
% positive scalar indicating a zoom factor
% negative scalar for default initial data
% 0 meaning to quit
%
% See also GO2, LO3
colors = ['b','g','r','c','m','y','k'];
ncolors = size(colors,2);
ic = 1;
global Ao bo
Ao = A; bo = [0;0;0];
tso = 25;
switch nargin
case 1, tsp = [0 tso]; tsm = -tsp;
case 2, tsp = [0 dt]; tsm = -tsp;
otherwise tsp = [0 dt]; tsm = -[0 et];
end
[E,D] = eig(A);
disp('Eigenvalues')
disp(diag(D))
disp(' ')
disp('Eigenvectors')
disp(E)
ax = 10 * [-1 1 -1 1 -1 1];
clf;
axis(ax);
zoom on;
rotate3d on;
stopit = 0;
while stopit == 0
y0 = input('Initial Data : ');
if y0 < 0,
y0 = [1 1 1; 1 1 -1; -1 1 1; -1 1 -1];
end
if size(y0,2) == 1
if y0 == 0
stopit = 1;
else
zoom(y0);
end
else
for i=1:size(y0,1)
[t,y] = ode45('lo3',tsp,y0(i,:));
hold on; plot3(y(:,1),y(:,2),y(:,3),colors(ic)); hold off;
[t,y] = ode45('lo3', tsm,y0(i,:));
hold on; plot3(y(:,1),y(:,2),y(:,3),colors(ic)); hold off;
ic = ic + 1; if ic > ncolors, ic = 1; end
end
end
end