JuliaDynamics / DynamicalSystems.jl

Award winning software library for nonlinear dynamics and nonlinear timeseries analysis
https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/
Other
849 stars 95 forks source link

Improving the logo #110

Open freemin7 opened 4 years ago

freemin7 commented 4 years ago

The current logo is acyclic. I think that using small, invisible nudges along the trajectory it might possible to move it into to a cyclic orbit.

If you find that interesting, could you provide me with the code that a generated the original animation?

Datseris commented 4 years ago

I also thought about it, and I thought of a more "natural" way to make it cyclic. We can evolve the pendulum until we reach a state that in the state space has very small distance from the initial condition: use integrator interface and do

integ = integrator(double_pendulum(u0))
while norm(integ.u - u0) > ε
step!(integ)
end

stopping the loop at that state and starting over might not have much visual impact. Of course, one will have to midfy the plotting code heavily. my only fear is that this might lead to huge evolution times integ.t. But on the other hand, we can choose a u0 that leads to very small ones.

But I also like a lot your idea of small nudges.

I'll try to find the plotting code and report back. The initial logo was actually done in PyPlot, because back then Julia didn't have any worthy plotting packages. Might be worthy to use MakieLayout now.

freemin7 commented 4 years ago

If you do not have much attachment to he current trajectory that would make searching a solution easier.

I have not solved non-linear trajectory optimization problems like this before. However i've been meaning to dig into these for a some time.

To increase the chances of finding a solution i was thinking of solving the equation forward and back backward in time. Since finding a trajectory under minimal control that leads upward is hard, while there needs to be some trajectory that got it into that state before which satisfies all the physical conditions without the presence of control.

freemin7 commented 4 years ago

If you have a certain "motive" of the double pendulum you really like, one could do multiple shooting and include those and their evolution, provided the total system energy is similar.

Datseris commented 4 years ago

If you do not have much attachment to he current trajectory that would make searching a solution easier.

Nope, not at all. In fact, I didn't even play around with the initial condition, just saved the first one.... :D

If you have a certain "motive" of the double pendulum you really like, one could do multiple shooting and include those and their evolution, provided the total system energy is similar

Something that is crucial is that at some point in the time evolution the disks align to form the Julia logo, like it happens in the static JuliaDynamics logo.

Datseris commented 3 years ago

So I now have a makie version of the double pendulum that looks like this:

https://user-images.githubusercontent.com/19669089/129259342-cd866d2d-b3e9-46b0-8429-a21be5a231ee.mp4

Code:

using DynamicalSystems
using DataStructures: CircularBuffer

# %% Simulate
const L1 = 1.0
const L2 = 0.9

M = 2
u0 = [π/3, 0, 3π/4, -2]
u0s = [u0 .+ [0, 0, 0.005*i, 0] for i in 0:M-1]
tail = 200
dp = Systems.double_pendulum(u0; L1, L2)

datas = trajectory.(Ref(dp), 200, u0s; Δt = 0.005)

function xycoords(data)
    θ1 = data[:, 1]
    θ2 = data[:, 3]
    x1 = L1 * sin.(θ1)
    y1 = -L1 * cos.(θ1)
    x2 = x1 + L2 * sin.(θ2)
    y2 = y1 - L2 * cos.(θ2)
    return x1,x2,y1,y2
end

# x1,x2,y1,y2 = xycoords(data)
coords = xycoords.(datas)

# %% Plot
fig = Figure(); display(fig)
ax = fig[1,1] = Axis(fig)
ax.aspect = DataAspect()
l = 1.05(L1+L2)
xlims!(ax, -l, l)
ylims!(ax, -l, 0.5l)
colors = reverse(COLORS[1:M])
hidespines!(ax)
hidedecorations!(ax)

function maketraj!(ax, x2, y2, c)
    traj = CircularBuffer{Point2f0}(tail)
    fill!(traj, Point2f0(x2, y2))
    traj = Observable(traj)
    tailcol =  [RGBAf0(c.r, c.g, c.b, (i/tail)^2) for i in 1:tail]
    lines!(ax, traj; linewidth = 3, color = tailcol)
    return traj
end
function makerod!(ax, c)
    rod = Observable([Point2f0(0, 0), Point2f0(0, -L1), Point2f0(0, -L1-L2)])
    node = Observable(Point2f0(0, -L1))
    ball = Observable(Point2f0(0, -L1-L2))
    lines!(ax, rod; linewidth = 4, color = lighten_color(c))
    scatter!(ax, node; marker = :circle, strokewidth = 2, strokecolor = lighten_color(c),
        color = darken_color(c, 2), markersize = 8
    )
    scatter!(ax, ball; marker = :circle, strokewidth = 3, strokecolor = lighten_color(c),
        color = darken_color(c, 2), markersize = 12
    )
    return rod, node, ball
end

trajs = [maketraj!(ax, coords[j][2][1], coords[j][4][1], colors[j]) for j in 1:M]
rods = [makerod!(ax, c) for c in colors]

function update!(trajs, rods, coords, i = 1)
    for j in 1:length(u0s)
        traj = trajs[j]
        rod, node, ball = rods[j]
        x1, x2, y1, y2 = coords[j]
        push!(traj[], Point2f0(x2[i], y2[i]))
        traj[] = traj[]
        rod[] = [Point2f0(0, 0), Point2f0(x1[i], y1[i]), Point2f0(x2[i], y2[i])]
        node[] = Point2f0(x1[i], y1[i])
        ball[] = Point2f0(x2[i], y2[i])
    end
end

update!(trajs, rods, coords)

# %% Animate
file = projectdir("video", "doublependulum.mp4")
mkpath(dirname(file))
record(fig, file; framerate = 60, compression = 10) do io
    for i = 1:3:4000
        update!(trajs, rods, coords, i)
        recordframe!(io)  # record a new frame
    end
end