Research

% Three anisotropic Gaussians, the building block of Gaussian splatting.
[X, Y] = meshgrid(linspace(-3, 3, 160));
G = zeros(size(X));
mu    = [-1.2 0.4; 0.8 -0.6; 1.3 1.1];
theta = [0.5 -0.9 1.6];
s     = [0.9 0.3; 0.7 0.25; 0.5 0.4];        % standard deviations along/across
w     = [1 0.8 0.6];

for k = 1:3
  c = cos(theta(k)); sn = sin(theta(k));
  dx = X - mu(k, 1); dy = Y - mu(k, 2);
  u = ( c*dx + sn*dy) / s(k, 1);
  v = (-sn*dx + c*dy) / s(k, 2);
  G = G + w(k) * exp(-0.5 * (u.^2 + v.^2));
end

contourf(X, Y, G, 14)
colormap magma
axis equal
title('Σ wₖ · exp(−½ (x−μₖ)ᵀ Σₖ⁻¹ (x−μₖ))')
xlabel('x'), ylabel('y')

sigma = 10; rho = 28; beta = 8/3;
f = @(tt, u) [sigma*(u(2) - u(1)); u(1)*(rho - u(3)) - u(2); u(1)*u(2) - beta*u(3)];
[tt, u] = ode45(f, [0 30], [1; 1; 1]);

plot3(u(:, 1), u(:, 2), u(:, 3), 'LineWidth', 0.7)
view(-25, 15)
title('Lorenz attractor')
subtitle(sprintf('σ = %g, ρ = %g, β = %.3g  ·  %d steps', sigma, rho, beta, numel(tt)))
xlabel('x'), ylabel('y'), zlabel('z')

% x ← r x (1 − x): period doubling into chaos
r = linspace(2.8, 4, 700);
x = 0.5 * ones(size(r));
for k = 1:300, x = r .* x .* (1 - x); end       % let transients die out
R = []; X = [];
for k = 1:80
  x = r .* x .* (1 - x);
  R = [R r]; X = [X x];
end
scatter(R, X, 1, 'k', 'filled')
xlabel('r'), ylabel('x')
title('Bifurcation diagram of the logistic map')

% Rosenbrock's banana valley and a gradient-descent path
f = @(x, y) (1 - x).^2 + 100*(y - x.^2).^2;
[X, Y] = meshgrid(linspace(-1.5, 1.5, 200), linspace(-0.5, 1.6, 200));
contour(X, Y, log10(f(X, Y) + 1e-3), 25)
colormap viridis
hold on
p = [-1.2 1]; path = p;
for k = 1:4000
  g = [-2*(1 - p(1)) - 400*p(1)*(p(2) - p(1)^2), 200*(p(2) - p(1)^2)];
  p = p - 0.0012 * g;
  if mod(k, 20) == 0, path = [path; p]; end
end
plot(path(:, 1), path(:, 2), 'r.-')
scatter(1, 1, 80, 'k', 'filled')
hold off
title(sprintf('after 4000 steps: (%.3f, %.3f)', p(1), p(2)))