I would like to simulate a Lorenz system, showing how the solution vector (x, y, z) changes with time. I know how to plot x, y and z data using PyPlot.plot3D, but how does one do a time plot, that is an animation based on t, x, y and z data one has? To be clear, I know how to solve the Lorenz equations in Julia, that part is easy, I just want to know how to create an animation with it in Julia.
Here is a video of an animation that is roughly along the lines of what I am asking how to do.
Thanks mate, the problem now is that this won’t seem to work. Here is my Julia script:
# This is largely copied from the ODE.jl repository
# Import package manager
using Pkg;
# Install and import ODE
Pkg.add("ODE")
using ODE;
function f(t, r)
# Extract the coordinates from the r vector
(x, y, z) = r
# The Lorenz equations
dx_dt = sigma*(y - x)
dy_dt = x*(rho - z) - y
dz_dt = x*y - bet*z
# Return the derivatives as a vector
[dx_dt; dy_dt; dz_dt]
end;
# Define time vector and interval grid
const dt = 0.001
const tf = 2e2
t = 0:dt:tf
# Initial position in space
const r0 = [-2.0; 3.0; 0.5]
# Constants sigma, rho and beta
const sigma = 10.0
const rho = 28.0
const bet = 8.0/3.0;
(t, pos) = ode78(f, r0, t)
x = map(v -> v[1], pos)
y = map(v -> v[2], pos)
z = map(v -> v[3], pos);
# Get PyPlot and load it
Pkg.add("PyPlot")
using PyPlot;
Pkg.add("Plots")
using Plots
PyPlot.figure(1)
PyPlot.plot3D(x, y, z);
PyPlot.xlabel("x")
PyPlot.ylabel("y")
PyPlot.zlabel("z")
PyPlot.figure(2)
PyPlot.plot3D(t, x, y);
PyPlot.xlabel("t")
PyPlot.ylabel("x")
PyPlot.zlabel("y")
PyPlot.figure(3)
PyPlot.plot3D(t, y, z);
PyPlot.xlabel("t")
PyPlot.ylabel("y")
PyPlot.zlabel("z")
PyPlot.figure(4)
PyPlot.plot3D(t, x, z);
PyPlot.xlabel("t")
PyPlot.ylabel("x")
PyPlot.zlabel("z")
fig, ax = subplots(1, 3, sharex=true, sharey=true, figsize=(16,8))
ax[1][:plot](x, y)
ax[1][:set_title]("X-Y cut")
ax[2][:plot](x, z)
ax[2][:set_title]("X-Z cut")
ax[3][:plot](y, z)
ax[3][:set_title]("Y-Z cut");
anim = @animate for i=1:length(t)
plot3D(x[i],y[i],z[i])
end
gif(anim, "/tmp/anim_fps15.gif", fps = 15)
Yet, I keep getting this error:
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.
│ caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.
│ caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.
│ caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.
│ caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.
│ caller = top-level scope at none:0
└ @ Core none:0
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.
│ caller = top-level scope at none:0
└ @ Core none:0
PyError ($(Expr(:escape, :(ccall(#= /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'TypeError'>
TypeError("object of type 'float' has no len()")
File "/usr/lib/python3.7/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 1530, in plot
zs = np.broadcast_to(zs, len(xs))
Stacktrace:
[1] pyerr_check at /home/fusion809/.julia/packages/PyCall/ttONZ/src/exception.jl:60 [inlined]
[2] pyerr_check at /home/fusion809/.julia/packages/PyCall/ttONZ/src/exception.jl:64 [inlined]
[3] macro expansion at /home/fusion809/.julia/packages/PyCall/ttONZ/src/exception.jl:84 [inlined]
[4] __pycall!(::PyCall.PyObject, ::Ptr{PyCall.PyObject_struct}, ::PyCall.PyObject, ::Ptr{Nothing}) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44
[5] _pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{Float64,Float64,Float64}, ::Int64, ::Ptr{Nothing}) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:29
[6] _pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{Float64,Float64,Float64}, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:11
[7] #pycall#110(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::PyCall.PyObject, ::Type{PyCall.PyAny}, ::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:86
[8] pycall(::PyCall.PyObject, ::Type{PyCall.PyAny}, ::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:86
[9] #plot3D#172(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyPlot/cdCMF/src/plot3d.jl:50
[10] plot3D(::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyPlot/cdCMF/src/plot3d.jl:48
[11] macro expansion at /data/GitHub/mine/math/julia-scripts/Lorenz.jl:82 [inlined]
[12] top-level scope at /home/fusion809/.julia/packages/Plots/oiirH/src/animation.jl:163
[13] include at ./boot.jl:326 [inlined]
[14] include_relative(::Module, ::String) at ./loading.jl:1038
[15] include(::Module, ::String) at ./sysimg.jl:29
[16] include(::String) at ./client.jl:403
[17] top-level scope at In[3]:1
To try and correct the error I was getting I added:
Pkg.add("Iterators")
using Iterators;
before the animation lines, but that just gave:
Unsatisfiable requirements detected for package Iterators [a4bce56a]:
Iterators [a4bce56a] log:
├─possible versions are: [0.1.0-0.1.10, 0.2.0, 0.3.0-0.3.1] or uninstalled
├─restricted to versions * by an explicit requirement, leaving only versions [0.1.0-0.1.10, 0.2.0, 0.3.0-0.3.1]
└─restricted by julia compatibility requirements to versions: uninstalled — no versions left
Stacktrace:
[1] #propagate_constraints!#61(::Bool, ::Function, ::Pkg.GraphType.Graph, ::Set{Int64}) at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/GraphType.jl:1007
[2] propagate_constraints! at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/GraphType.jl:948 [inlined] (repeats 2 times)
[3] #simplify_graph!#121(::Bool, ::Function, ::Pkg.GraphType.Graph, ::Set{Int64}) at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/GraphType.jl:1462
[4] simplify_graph! at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/GraphType.jl:1462 [inlined] (repeats 2 times)
[5] resolve_versions!(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}, ::Nothing) at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/Operations.jl:371
[6] resolve_versions! at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/Operations.jl:315 [inlined]
[7] #add_or_develop#63(::Array{Base.UUID,1}, ::Symbol, ::Function, ::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}) at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/Operations.jl:1172
[8] #add_or_develop at ./none:0 [inlined]
[9] #add_or_develop#17(::Symbol, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}) at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/API.jl:59
[10] #add_or_develop at ./none:0 [inlined]
[11] #add_or_develop#16 at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/API.jl:36 [inlined]
[12] #add_or_develop at ./none:0 [inlined]
[13] #add_or_develop#13 at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/API.jl:34 [inlined]
[14] #add_or_develop at ./none:0 [inlined]
[15] #add_or_develop#12(::Base.Iterators.Pairs{Symbol,Symbol,Tuple{Symbol},NamedTuple{(:mode,),Tuple{Symbol}}}, ::Function, ::String) at /build/julia/src/julia-1.1.1/usr/share/julia/stdlib/v1.1/Pkg/src/API.jl:33
[16] add(::String) at ./none:0
[17] top-level scope at none:0
[18] include at ./boot.jl:326 [inlined]
[19] include_relative(::Module, ::String) at ./loading.jl:1038
[20] include(::Module, ::String) at ./sysimg.jl:29
[21] include(::String) at ./client.jl:403
[22] top-level scope at In[4]:1
and now I am kind of stuck. Could you help perhaps?
I’ve also tried this script (keeping it simple, in essence):
# This is largely copied from the ODE.jl repository
# Import package manager
using Pkg;
# Install and import ODE
Pkg.add("ODE")
using ODE;
function f(t, r)
# Extract the coordinates from the r vector
(x, y, z) = r
# The Lorenz equations
dx_dt = sigma*(y - x)
dy_dt = x*(rho - z) - y
dz_dt = x*y - bet*z
# Return the derivatives as a vector
[dx_dt; dy_dt; dz_dt]
end;
# Define time vector and interval grid
const dt = 0.001
const tf = 2e2
t = 0:dt:tf
# Initial position in space
const r0 = [-2.0; 3.0; 0.5]
# Constants sigma, rho and beta
const sigma = 10.0
const rho = 28.0
const bet = 8.0/3.0;
(t, pos) = ode78(f, r0, t)
x = map(v -> v[1], pos)
y = map(v -> v[2], pos)
z = map(v -> v[3], pos);
# Get PyPlot and load it
Pkg.add("PyPlot")
using PyPlot;
#Pkg.add("Plots")
#Pkg.add("Iterators")
#using Plots
#PyPlot.figure(1)
#PyPlot.plot3D(x, y, z);
#PyPlot.xlabel("x")
#PyPlot.ylabel("y")
#PyPlot.zlabel("z")
#PyPlot.figure(2)
#PyPlot.plot3D(t, x, y);
#PyPlot.xlabel("t")
#PyPlot.ylabel("x")
#PyPlot.zlabel("y")
#PyPlot.figure(3)
#PyPlot.plot3D(t, y, z);
#PyPlot.xlabel("t")
#PyPlot.ylabel("y")
#PyPlot.zlabel("z")
#PyPlot.figure(4)
#PyPlot.plot3D(t, x, z);
#PyPlot.xlabel("t")
#PyPlot.ylabel("x")
#PyPlot.zlabel("z")
#fig, ax = subplots(1, 3, sharex=true, sharey=true, figsize=(16,8))
#ax[1][:plot](x, y)
#ax[1][:set_title]("X-Y cut")
#ax[2][:plot](x, z)
#ax[2][:set_title]("X-Z cut")
#ax[3][:plot](y, z)
#ax[3][:set_title]("Y-Z cut");
anim = @animate for i=1:length(t)
plot3D(x[i],y[i],z[i])
end
gif(anim, "/tmp/anim_fps15.gif", fps = 15)
it too failed with the error:
PyError ($(Expr(:escape, :(ccall(#= /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'TypeError'>
TypeError("object of type 'float' has no len()")
File "/usr/lib/python3.7/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 1530, in plot
zs = np.broadcast_to(zs, len(xs))
Stacktrace:
[1] pyerr_check at /home/fusion809/.julia/packages/PyCall/ttONZ/src/exception.jl:60 [inlined]
[2] pyerr_check at /home/fusion809/.julia/packages/PyCall/ttONZ/src/exception.jl:64 [inlined]
[3] macro expansion at /home/fusion809/.julia/packages/PyCall/ttONZ/src/exception.jl:84 [inlined]
[4] __pycall!(::PyCall.PyObject, ::Ptr{PyCall.PyObject_struct}, ::PyCall.PyObject, ::Ptr{Nothing}) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44
[5] _pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{Float64,Float64,Float64}, ::Int64, ::Ptr{Nothing}) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:29
[6] _pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{Float64,Float64,Float64}, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:11
[7] #pycall#110(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::PyCall.PyObject, ::Type{PyCall.PyAny}, ::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:86
[8] pycall(::PyCall.PyObject, ::Type{PyCall.PyAny}, ::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:86
[9] #plot3D#172(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyPlot/cdCMF/src/plot3d.jl:50
[10] plot3D(::Float64, ::Vararg{Float64,N} where N) at /home/fusion809/.julia/packages/PyPlot/cdCMF/src/plot3d.jl:48
[11] macro expansion at /data/GitHub/mine/math/julia-scripts/Lorenz.jl:83 [inlined]
[12] top-level scope at /home/fusion809/.julia/packages/Plots/oiirH/src/animation.jl:163
[13] include at ./boot.jl:326 [inlined]
[14] include_relative(::Module, ::String) at ./loading.jl:1038
[15] include(::Module, ::String) at ./sysimg.jl:29
[16] include(::String) at ./client.jl:403
[17] top-level scope at In[8]:1
I don’t see the connection between your code and the Plots documentation - note that Plots and PyPlot are separate plotting packages, and while PyPlot is one of the backends of Plots, you shouldn’t be calling both at the same time.
The Plots docs actually have the Lorenz attractor as an example on their [homepage[(Home · Plots):
# define the Lorenz attractor
mutable struct Lorenz
dt; σ; ρ; β; x; y; z
end
function step!(l::Lorenz)
dx = l.σ*(l.y - l.x) ; l.x += l.dt * dx
dy = l.x*(l.ρ - l.z) - l.y ; l.y += l.dt * dy
dz = l.x*l.y - l.β*l.z ; l.z += l.dt * dz
end
attractor = Lorenz((dt = 0.02, σ = 10., ρ = 28., β = 8//3, x = 1., y = 1., z = 1.)...)
# initialize a 3D plot with 1 empty series
plt = plot3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
title = "Lorenz Attractor", marker = 2)
# build an animated gif by pushing new points to the plot, saving every 10th frame
@gif for i=1:1500
step!(attractor)
push!(plt, attractor.x, attractor.y, attractor.z)
end every 10
``Using plots was an attempt to fix the error shown when I didn’t use it. I thought it might be necessary given that it was used later in the docs you showed me.
Thanks for the code!
A bit late to the party, but I am just starting to learn Julia and came across the Lorenz attractor example in the docs. While the approach is quite intuitive and elegant, I am struggling to extend this to subplots, i.e. two trajectories of the attractor next to each other.
I have tried to push to the subplots instead of the overall plot, but as far as I understand the error message, subplots do not support the pushing of new data?
N = 200
u1 = 20*rand(3, N)
u2 = 15*ones(3, N)
plt = plot3d(
1,
xlim = (-30, 30),
ylim = (-30, 30),
zlim = (0, 60),
title = ["Some trajectory" "Other trajectory"],
legend = false,
marker = 2,
layout = (1, 2)
)
@gif for i in range(1, size(u1, 2))
push!(plt[1], u1[:, i]...)
push!(plt[2], u2[:, i]...)
end every 10
Giving me the following error:
MethodError: no method matching push!(::Plots.Subplot{Plots.GRBackend}, ::Float64)
Closest candidates are:
push!(::Any, ::Any, !Matched::Any)
@ Base abstractarray.jl:3389
push!(::Any, ::Any, !Matched::Any, !Matched::Any...)
@ Base abstractarray.jl:3390
push!(!Matched::StatsBase.AbstractHistogram{T, 1}, ::Real) where T
@ StatsBase some_path\.julia\packages\StatsBase\WLz8A\src\hist.jl:294
...
Stacktrace:
[1] push!(A::Plots.Subplot{Plots.GRBackend}, a::Float64, b::Float64)
@ Base .\abstractarray.jl:3389
[2] push!(A::Plots.Subplot{Plots.GRBackend}, a::Float64, b::Float64, c::Float64)
@ Base .\abstractarray.jl:3390
[3] macro expansion
@ some_path\MMDS_Julia\ex2.ipynb:19 [inlined]
[4] top-level scope
@ some_path\.julia\packages\Plots\sxUvK\src\animation.jl:251
Any suggestion how to achieve the wanted effect?
Welcome to the forum @CM1 !
It’s probably slower since the combined plot has to be recreated every step of the animation, but this worked for me:
using Plots
N = 200
u1 = 20*rand(3, N)
u2 = 15*ones(3, N)
plt1 = plot3d(
1,
xlim = (-30, 30),
ylim = (-30, 30),
zlim = (0, 60),
title = ["Some trajectory"],
legend = false,
marker = 2,
)
plt2 = plot3d(
1,
xlim = (-30, 30),
ylim = (-30, 30),
zlim = (0, 60),
title = ["Other trajectory"],
legend = false,
marker = 2,
)
@gif for i in range(1, size(u1, 2))
push!(plt1, u1[:, i]...)
push!(plt2, u2[:, i]...)
plot(plt1, plt2)
end every 10
The docstring of plot mentions this about extracting subplots, but I couldn’t get it to work in the intended way (creating a new plot seems to remove the previously drawn points):
help?> plot
...
Extract a subplot from an existing plot.
Examples
≡≡≡≡≡≡≡≡≡≡
julia> p1, p2 = plot(1:2), plot(10:20)
julia> pl = plot(p1, p2) # plot containing 2 subplots
julia> plot(pl.subplots[1]) # extract 1st subplot as a standalone plot
julia> plot(pl.subplots[2]) # extract 2nd subplot as a standalone plot
Thank you! It seems like an overkill to me to create two distinct plot objects, but it’s working like a charm 
That’s great!
Yeah, would be nice to know if there is a way to modify subplots in place. I don’t know how (but I’m also not an expert with Plots.jl). Perhaps you can open a separate thread about it here 
EDIT: And I’m not sure how big the overhead actually is, since only the “wrapper” plot gets recreated at every step. The subplots are just modified.
You’re right, at this scale it doesn’t seem like a huge overhead, but I am not sure what happens under the hood, maybe someone can clarify whether this is just a cosmetic issue or actually slower.
I created a new discussion in this thread.