How to create an animation based on t, x, y and z data?

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.

1 Like

https://docs.juliaplots.org/latest/animations/

1 Like

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[(http://docs.juliaplots.org/latest/):

# 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
``
1 Like

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!