%% Lesson 10, extra -- Chaos!
% Based on work for PH429: Chaos & Nonlinear Dynamics. The following was
% translated from the original NumPy/Matplotlib code. See Strogatz' Nonlinear
% Dynamics and Chaos (the source for much of this information) if you want to
% know more!

close all; clear; clc;

% Simple systems of nonlinear differential equations (also iterated nonlinear
% maps) often exhibit surprisingly complex behavior. Lorenz's "strange
% attractor," discovered in 1963, was an example of this: he discovered that an
% exceedingly simplified (and entirely deterministic!) model of the weather was
% still able to produce unpredictable fluctuations in temperature, wind speed,
% etc. -- the system exponentially amplified miniscule parameter changes, so
% every decimal place mattered. Further work with analog computers and the
% then-emerging digital computers led to the field we now call chaos theory.

% Chaos has much to do with fractal geometry, so we'll start by generating a
% few fractals. Because... why not?

%% The Chaos Game & Sierpinski Triangles
% The game is simple:
% 1. Place three vertices in the plane, corresponding to the vertices of an
%    equilateral triangle.
% 2. Randomly choose a point P(0).
% 3. Each turn, choose a random vertex, and plot a point P(n) halfway between
%    that vertex and point P(n - 1).
% Let's see what happens.

n1 = 5e5; % points
vertices = [0, 0; .5, .5*sqrt(3); 1, 0];
r = .5; % fraction of the way between points and vertices (here halfway)

% point coordinates
points = zeros(n1, 2);
choices = randi([1 3], n1, 1);

% randomly choose the first point somewhere in the plane
points(1, :) = 100*randn(1, 2);

% choose the rest
% note: the loop here is actually necessary, like it was in G-S -- each
% iteration depends upon the previous one.
for i = 2:n1
	choice = choices(i);
	vertex = vertices(choice, :);
	last_point = points(i - 1, :);

	points(i, :) = last_point + r*(vertex - last_point);
end

% plot the thing
figure;
scatter(points(:, 1), points(:, 2), 1);
xlim([-.1 1.1]);
ylim([-.2 1]);
xticks([]);
yticks([]);
pbaspect([1 1 1]); % make the plot box square
title('The Sierpinski Triangle, by Chance', 'Interpreter', 'latex');

% A rather important property of fractals is their fractional dimension -- this
% triangle, for example, is not quite two-dimensional but is clearly more than
% one-dimensional. So how do we measure this dimension?

% One simple measure, known as self-similarity dimension, is relatively easy to
% find for the Sierpinski triangle. The idea is we simply examine how many
% smaller versions of a given shape could fit into that shape, and what size
% they would be in relation. For example, a square with scale 1 can be made up
% of 4 smaller squares scaled by 1/2, or 9 smaller squares with scale 1/3:
% -----------------     -----------------
% |               |     |       |       |
% |               |     |  1/2  |  1/2  |
% |               |     |       |       |
% |       1       |  =  |-------|-------|
% |               |     |       |       |
% |               |     |  1/2  |  1/2  |
% |               |     |       |       |
% -----------------     -----------------
% A square thus has self-similarity dimension log₂4 = log₃9 = 2 (as one might
% expect). The Sierpinski triangle is composed of 3 copies of itself scaled by
% 1/2, and thus has self-similarity dimension log₂3 ≈ 1.585.

% I wonder if we could have MATLAB calculate these...

%% Iterated Maps
% One way to create chaotic behavior is to iterate a nonlinear map: e.g., plot
% the set {x, sin(x), sin(sin(x)), sin(sin(sin(x))), ...}. The values will
% generally hop unpredictably around, but will often stay closer to some values
% than others...

% The following plot is a visual representation of these iterations, called a
% cobweb. Cobwebs are often used to qualitatively analyze nonlinear maps. The
% idea is, since the output of the map will be fed back into the input, we plot
% a point on the map's curve, translate horizontally to the line y = x, then
% translate vertically back to the map and repeat. This gives a visual
% representation of plugging the output back into the input.

% When you see this, think about which points on the map will "attract" the
% cobweb, which will "repel" it, and how one might go about making that
% distinction.

n2 = 10;
map = @(x) 2*sin(x);

% create the map curve and the center line
x = linspace(0, pi, 100).';
y_map = map(x);
y_line = x;

% create the cobweb
map_pts = zeros(2*n2 + 1, 2);
map_pts(1, :) = [.2, 0];
for i = 1:n2
	last_pt = map_pts(2*i - 1, :);
	mapped = [last_pt(1), map(last_pt(1))];
	on_the_line = [mapped(2), mapped(2)];

	map_pts(2*i, :) = mapped;
	map_pts(2*i + 1, :) = on_the_line;
end

% plot it
figure;
hold on;
plot(x, y_map, 'b');
plot(x, y_line, 'k');
plot(map_pts(:, 1), map_pts(:, 2), 'r');

start = map_pts(1, :);
scatter(start(1), start(2), 'r');
text(start(1) + .03, start(2) + .03, 'Start Point');

xlim([0 pi]);
xticks([0 pi/2 pi]);
ylim([0 2.1]);
xticklabels({'0', '\pi/2', '\pi'});

title('Cobweb of \(2 \sin x\)', 'Interpreter', 'latex');

%% Orbit Diagrams
% The above example was... not too crazy. The map x → α sin x *can* exhibit
% more complex behavior, but it requires α to be greater than π (so that we can
% escape the [0, π] interval). The map will first stay close to two points,
% then a few, then... altogether too many. This is usually plotted with what is
% called an orbit diagram, showing roughly how much time the map spends at
% certain values for a given α.

map_alpha = @(alpha, x) alpha.*sin(x);

n3 = 6e2; % iterations to go through
frames = 9e2; % run this many parallel alpha computations
toss = 2e2; % delete this many points from the start

min_alpha = 2.2;
max_alpha = 3.1;
alphas = linspace(min_alpha, max_alpha, frames);

% at least we can run the frames in parallel...
values = zeros(n3, frames);
values(1, :) = pi/4; % the first input each time
for i = 2:n3
	values(i, :) = map_alpha(alphas, values(i - 1, :));
end

% plot this crazy thing
values = values((toss + 1):end, :);
figure;
scatter(alphas.', values, 1, [.7 0 0], 'MarkerEdgeAlpha', .3);

title('Orbit of \(\alpha \sin x\)', 'Interpreter', 'latex');
xlabel('\alpha');
ylabel('map value');

%% The Phase Plane
% In continuous time, we have differential equations rather than iterated maps.
% One way to analyze higher-order equations is using the so-called phase plane,
% which is a plot of our state variable and its derivative. The idea is that a
% second-order differential equation really acts like a system in two
% dimensions, and is best analyzed directly as such.

% We'll analyze the motion of an damped pendulum, which has the governing
% equation d²θ/dt² + b dθ/dt + (g/L) sin θ = 0. For simplicity (it won't change
% the qualitative behavior), we'll let g/L = 1. Then (the phase plane bit)
% we break this into two differential equations,
%  * dθ/dt = ω and
%  * dω/dt = -(sin θ + bω),
% and make a plot of the *vector* (dθ/dt, dω/dt) in θ-ω space. We can easily
% plot the vector field with `quiver`, and MATLAB's `streamline` is rather
% effective at plotting either single trajectories or a whole suite of possible
% solutions.

b = .2;
min_theta = -3*pi;
max_theta = 3*pi;
min_omega = -3*pi;
max_omega = 3*pi;
theta_range = linspace(min_theta, max_theta, 300);
omega_range = linspace(min_omega, max_omega, 200);
[theta, omega] = meshgrid(theta_range, omega_range);
theta_dot = omega;
omega_dot = -(sin(theta) + b*omega);

% to cut down on the # of points we plot
d = 10;
[sparse_x, sparse_y] = ...
	meshgrid(d:d:length(theta_range), d:d:length(omega_range));
sparse = sub2ind(size(theta), sparse_y(:), sparse_x(:));

% plot this thing
figure;
hold on;
quiver(theta(sparse), omega(sparse), ...
	theta_dot(sparse), omega_dot(sparse), 'm');
streamline(theta, omega, theta_dot, omega_dot, ...
	theta(sparse_y, [1 end]), omega(sparse_y, [1 end]));

xlim([min_theta max_theta]);
xticks(linspace(min_theta, max_theta, 7));
ylim([min_omega max_omega]);

xlabel('pendulum angle');
ylabel('angular velocity');
title('Behavior of a Damped Pendulum', 'Interpreter', 'latex');

%% Strange Attractors
% I won't go into detail here -- suffice it to say the behavior of this system
% would be quite opaque without numerical integration. Fortunately, we have
% just that! Lorenz managed to prove analytically that this system had no
% closed cycles and yet all its trajectories eventually stayed within a bounded
% region of space, but beyond that... well, we'll see.
sigma = 10;
r = 28;
b = 8/3;

x_dot = @(x, y, z) sigma*(y - x);
y_dot = @(x, y, z) r*x - y - x.*z;
z_dot = @(x, y, z) x.*y - b*z;

[x, y, z] = meshgrid(linspace(-50, 50, 100));
u = x_dot(x, y, z);
v = y_dot(x, y, z);
w = z_dot(x, y, z);

figure;
hold on;
streamline(x, y, z, u, v, w, 0, 1, 0, [.1 2e4]);
scatter3(0, 1, 0);
view(3);
grid on;

text(0, 0, -2, 'Start Point');
xlabel('x');
ylabel('y');
zlabel('z');
title("Lorenz's Strange Attractor", 'Interpreter', 'latex');