Consider the well-known chaotic model of Henon & Heiles (motivated by the motion

Question

Consider the well-known chaotic model of Henon & Heiles (motivated by the motion of stars in a galaxy). The model is dp_1/dt = -q1 - 2q- q_1^2 + q_2^2, p_2(0) = 0.24, dq_1/dt = p_1, q_1(0) = 0.24, dq_2/dt = p_2, q_2(0) = 0.24, Integrate numerically this initial value problem using Matlab's ode113 solver with relative tolerance equal to 100 times eps and absolute tolerance equal to 10^-14. Integrate on the time Interval [0, 100]. Plot on the same graph the evolution in time of p_1, p_2, q_1 and q_2, respectively. Integrate on the time interval [0, 2 times 10^4]. Create the Poincare map of the system by taking q_1(t) = 0 and plotting q_2(t) against p_2(t). Only plot q_2(t) against p_2(t) at the points where q_1(t) = 0 and q_1(t) is increasing. The plot should use discrete values (dots) rather than the interpolated output (lines).

Explanation / Answer

File code.m

clear
clc

%% Part (a)

options = odeset('RelTol', 100*eps, 'AbsTol', 1e-14);
[t,y] = ode113(@derivatives, [0 100], [0.24; 0.24; 0.24; 0.24], options);

figure
plot(y(:,1)); hold on
plot(y(:,2)); hold on
plot(y(:,3)); hold on
plot(y(:,4));
legend('p_{1}', 'p_{2}', 'q_{1}', 'q_{2}');
axis([0 1000 -1 1]); xlabel('t');

%% Part (b)

options = odeset(options, 'Events', @q1Events);
[t,y,te,ye,ie] = ode113(@derivatives, [0 2*(10^4)], [0.24; 0.24; 0.24; 0.24], options);
figure, scatter(ye(:,2), ye(:,4));
xlabel('p_{2}'); ylabel('q_{2}');
title('q_{2} vs. p_{2}');

File derivatives.m

function dydt = derivatives(t, y)
    dydt = [-y(3)-2*y(3)*y(4);
            -y(4)-y(3)*y(3)+y(4)*y(4);
            y(1);
            y(2)];
end

File q1Events.m

function [position,isterminal,direction] = q1Events(t, y)
    position = y(3);
    isterminal = 0;
    direction = 1;
end

Keep all three files in the same folder and run the file code.m

Related Questions

Navigate

• Browse (All)
• Browse C
• Subjects

---