Animates the evolution of an initially tight group of points ... my intro to Julia

lorenz-animator.jl

using DifferentialEquations
using Plots
using Statistics
using LinearAlgebra

function lorenz!(du, u, p, t)
    x, y, z = u
    σ, ρ, β = p

    du[1] = σ * (y - x)
    du[2] = x * (ρ - z) - y
    du[3] = x * y - β * z
end

function track()
    # start near a known point
    u0 = [-7, 7, 25] + randn(3) * 1e-4
    tspan = (0.0, 25)
    prob = ODEProblem(lorenz!, u0, tspan, parammap = (10, 28, 8 / 3))
    solve(prob)
end

function center(x)
    x .- mean(x)
end

"""
Keeps track of the scale that a variable has shown
"""
function filter!(u, state, horizon=200)
    old, w = state
    z = reduce(max, abs.(u))
    old = max(old, z)
    z = (z + w*old) / (w + 1)
    state[1] = z
    state[2] = min(horizon, w + 1)
    1.1 * z
end

"""
Compute nice ticks for a graph that grows and shrinks.
"""
function ticks(xmin, xmax=abs(xmin))
    if xmin >= 0
        xmin = -xmax
    end
    biggest = max(abs(xmin), abs(xmax))
    scale = 10^ceil(log10(biggest)) .* [0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5]
    base = [tick for tick in scale if (xmin ≤ tick ≤ xmax) && (xmax / tick < 5)]
    vcat(reverse(-base), 0, base)
end

# Solve 20 different evolutions
balls = [track() for i in 1:20]
# turn off local rendering of graphs
ENV["GKSwstype"] = "nul"
# initial tracker state has no information
tx1 = [0.0, 0.0]
# we compute a projection of the attractor that lays things out well
# but is reasonable right side up
V = svd(transpose(hcat(balls[1].(LinRange(0,10,25))...))).V
up = transpose(V[:, 1:2]) * [0, 1, 0]
project = transpose(V[:, 1:2] * [up .* [-1, 1] reverse(up)])

# animate the solutions we got above
anim = @gif for t in LinRange(0, 15,1500)
    l = @layout [ a{0.40w} b{0.5w}  ]
    x = hcat([b(t) for b in balls]...)
    # plots the 3d projection
    p1 = plot(balls[1], vars = (1,2,3), legend = false, xlim = (-22,22), ylim = (-22,22), zlim = (0,45))
    # with added balls for the solutions
    scatter!(p1, x[1, :], x[2, :], x[3, :])
    
    # reset to center of mass 
    diffs = transpose([center(x[1, :]) center(x[2, :]) center(x[3, :])])
    projected = project * diffs
    px, py = projected[1, :], projected[2, :]
    xmax = filter!([px; py], tx1)
    # plut the solution after projection
    p2 = scatter(px, py, xlim = (-xmax, xmax), ylim = (-xmax, xmax), legend = false, ticks = ticks(-xmax,xmax))
    
    # lay down the animation frame
    plot(p1, p2, layout = l)
end