Newer
Older
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
% main
% Initial conditions
t0 = 0;
w0 = [1,1,1];
% Maximum time of integration
T = 100;
% Step definitions
h0 = 0.01; % initial step size
h = h0; % used step size
hmin = (1e-5)*h0; % minimum step size allowed
% Model parameters and definitions
sigma = 10;
rho = 28;
beta = 8/3;
F = @(s,w) Lorenz(s,w,sigma,rho,beta);
% Poincare Section definition
Section = @ (w) w(3) - 15;
tsec = [];
wsec = [];
% Result allocation
t = zeros(100000,1);
w = zeros(100000,max(size(w0)));
t(1) = t0;
w(1,:) = w0;
counter = 2;
% Initiate figure object
figure;
az = 0; el = 0; % viewpoint definition
%%%%%% ITERATE
while t(counter) <= T
% iterate a step
[t(counter), w(counter,:)] = RK4(t(counter-1),w(counter-1,:),h,F);
% evaluete if solutuion passes through the poincare section
if or(sign(Section(w(counter,:))) > sign(Section(w(counter-1,:))), sign(Section(w(counter,:))) < sign(Section(w(counter-1,:))))
h = h/2;
if h <= hmin
h = h0;
tsec = [tsec;t(counter)];
wsec = [wsec;w(counter,:)];
counter = counter + 1;
end
else
% plot result
plot3(w(1:counter,1),w(1:counter,2),w(1:counter,3))
xlabel('x'); ylabel('y'); zlabel('z')
hold on;
% plot points that pass through Poincare section
if max(size(wsec))>0
plot3(wsec(:,1),wsec(:,2),wsec(:,3),'o');
end
% plot current point in trajectory
plot3(w(counter,1),w(counter,2),w(counter,3),'*')
hold off
% rotate and view solution
az = mod(az + 0.5,360);
view(az,el)
pause(0.005);
%iterate
counter = counter + 1;
end
end
% eliminate unused entries in solution vector
t(counter+1:end) = [];
w(counter+1,:) = [];