Poincaré Integral Invariant Example
2nd Poincaré Invariant
Load modules
using GeometricIntegrators
using PoincareInvariants
using PyPlot[ Info: Installing matplotlib via the Conda matplotlib package...
[ Info: Running `conda install -q -y matplotlib` in root environment
Collecting package metadata (current_repodata.json): ...working... done
Solving environment: ...working... done
## Package Plan ##
environment location: /home/runner/.julia/conda/3
added / updated specs:
- matplotlib
The following packages will be downloaded:
package | build
---------------------------|-----------------
brotli-1.0.9 | he6710b0_2 375 KB
cycler-0.10.0 | py39h06a4308_0 16 KB
dbus-1.13.18 | hb2f20db_0 504 KB
expat-2.4.1 | h2531618_2 168 KB
fontconfig-2.13.1 | h6c09931_0 250 KB
fonttools-4.25.0 | pyhd3eb1b0_0 632 KB
freetype-2.10.4 | h5ab3b9f_0 596 KB
glib-2.69.1 | h5202010_0 1.7 MB
gst-plugins-base-1.14.0 | h8213a91_2 4.9 MB
gstreamer-1.14.0 | h28cd5cc_2 3.2 MB
icu-58.2 | he6710b0_3 10.5 MB
jpeg-9b | h024ee3a_2 214 KB
kiwisolver-1.3.1 | py39h2531618_0 80 KB
lcms2-2.12 | h3be6417_0 312 KB
libpng-1.6.37 | hbc83047_0 278 KB
libtiff-4.2.0 | h85742a9_0 502 KB
libuuid-1.0.3 | h1bed415_2 15 KB
libwebp-base-1.2.0 | h27cfd23_0 437 KB
libxcb-1.14 | h7b6447c_0 505 KB
libxml2-2.9.12 | h03d6c58_0 1.2 MB
lz4-c-1.9.3 | h295c915_1 185 KB
matplotlib-3.4.2 | py39h06a4308_0 26 KB
matplotlib-base-3.4.2 | py39hab158f2_0 5.6 MB
munkres-1.1.4 | py_0 13 KB
olefile-0.46 | py_0 33 KB
openjpeg-2.4.0 | h3ad879b_0 331 KB
pcre-8.45 | h295c915_0 207 KB
pillow-8.3.1 | py39h2c7a002_0 648 KB
pyparsing-2.4.7 | pyhd3eb1b0_0 59 KB
pyqt-5.9.2 | py39h2531618_6 4.7 MB
python-dateutil-2.8.2 | pyhd3eb1b0_0 233 KB
qt-5.9.7 | h5867ecd_1 68.5 MB
sip-4.19.13 | py39h2531618_0 279 KB
tornado-6.1 | py39h27cfd23_0 595 KB
zstd-1.4.9 | haebb681_0 480 KB
------------------------------------------------------------
Total: 108.1 MB
The following NEW packages will be INSTALLED:
brotli pkgs/main/linux-64::brotli-1.0.9-he6710b0_2
cycler pkgs/main/linux-64::cycler-0.10.0-py39h06a4308_0
dbus pkgs/main/linux-64::dbus-1.13.18-hb2f20db_0
expat pkgs/main/linux-64::expat-2.4.1-h2531618_2
fontconfig pkgs/main/linux-64::fontconfig-2.13.1-h6c09931_0
fonttools pkgs/main/noarch::fonttools-4.25.0-pyhd3eb1b0_0
freetype pkgs/main/linux-64::freetype-2.10.4-h5ab3b9f_0
glib pkgs/main/linux-64::glib-2.69.1-h5202010_0
gst-plugins-base pkgs/main/linux-64::gst-plugins-base-1.14.0-h8213a91_2
gstreamer pkgs/main/linux-64::gstreamer-1.14.0-h28cd5cc_2
icu pkgs/main/linux-64::icu-58.2-he6710b0_3
jpeg pkgs/main/linux-64::jpeg-9b-h024ee3a_2
kiwisolver pkgs/main/linux-64::kiwisolver-1.3.1-py39h2531618_0
lcms2 pkgs/main/linux-64::lcms2-2.12-h3be6417_0
libpng pkgs/main/linux-64::libpng-1.6.37-hbc83047_0
libtiff pkgs/main/linux-64::libtiff-4.2.0-h85742a9_0
libuuid pkgs/main/linux-64::libuuid-1.0.3-h1bed415_2
libwebp-base pkgs/main/linux-64::libwebp-base-1.2.0-h27cfd23_0
libxcb pkgs/main/linux-64::libxcb-1.14-h7b6447c_0
libxml2 pkgs/main/linux-64::libxml2-2.9.12-h03d6c58_0
lz4-c pkgs/main/linux-64::lz4-c-1.9.3-h295c915_1
matplotlib pkgs/main/linux-64::matplotlib-3.4.2-py39h06a4308_0
matplotlib-base pkgs/main/linux-64::matplotlib-base-3.4.2-py39hab158f2_0
munkres pkgs/main/noarch::munkres-1.1.4-py_0
olefile pkgs/main/noarch::olefile-0.46-py_0
openjpeg pkgs/main/linux-64::openjpeg-2.4.0-h3ad879b_0
pcre pkgs/main/linux-64::pcre-8.45-h295c915_0
pillow pkgs/main/linux-64::pillow-8.3.1-py39h2c7a002_0
pyparsing pkgs/main/noarch::pyparsing-2.4.7-pyhd3eb1b0_0
pyqt pkgs/main/linux-64::pyqt-5.9.2-py39h2531618_6
python-dateutil pkgs/main/noarch::python-dateutil-2.8.2-pyhd3eb1b0_0
qt pkgs/main/linux-64::qt-5.9.7-h5867ecd_1
sip pkgs/main/linux-64::sip-4.19.13-py39h2531618_0
tornado pkgs/main/linux-64::tornado-6.1-py39h27cfd23_0
zstd pkgs/main/linux-64::zstd-1.4.9-haebb681_0
Preparing transaction: ...working... done
Verifying transaction: ...working... done
Executing transaction: ...working... doneDefine functions for phasespace surface, one- and two-form and vector field
function surface(s, t, offset=1.0, factor=1.0)
[(offset+s)/factor, (offset+t)/factor, s, t]
end
function theta(t, q)
return [x[3] - x[2], x[4], 0., 0.]
end
function omega(t, q, B)
B .= [[ 0. 1. -1. 0.]
[-1. 0. 0. -1.]
[ 1. 0. 0. 0.]
[ 0. 1. 0. 0.]]
end
function f(t, q, f)
f[1] = q[3]
f[2] = q[4]
f[3] = +q[4]
f[4] = -q[3]
endSet parameters
const nd = 4
const n₀ = 100
const nt = 200
const Δt = 0.1Trapezoidal Quadrature
Create Poincaré Invariant
pinv = PoincareInvariant2ndTrapezoidal(
x₀ -> ode = ODE(f, x₀), surface, omega,
Δt, nd, n₀, n₀, nt)Compute solution
tab = TableauGLRK(1)
int = Integrator(pinv.equ, tab, pinv.Δt)
sol = integrate(pinv.equ, int, pinv.ntime)
nothing # hidePlot solution
fig = figure(figsize=(4,4))
plot(sol.q[1,0,:], sol.q[2,0,:], ".", markersize=.5)
savefig("example_2nd_trapezoidal_initial_state.svg")
fig = figure(figsize=(4,4))
plot(sol.q[1,end,:], sol.q[2,end,:], ".", markersize=.5)
savefig("example_2nd_trapezoidal_final_state.svg")
fig = figure(figsize=(4,4))
for i in 1:n₀
plot3D(sol.q[1,:,i], sol.q[2,:,i], collect(0:nt)*Δt)
end
savefig("example_2nd_trapezoidal_evolution1.png")
fig = figure(figsize=(4,4))
for i in 0:20:nt
plot3D(sol.q[1,i,:], sol.q[2,i,:], i*Δt)
end
savefig("example_2nd_trapezoidal_evolution2.png")
Compute Poincarè invariant
I, J, ΔI, ΔJ = evaluate_poincare_invariant(pinv, sol)
pinv.I[0]Plot invariant error
yf = matplotlib[:ticker][:ScalarFormatter]()
yf[:set_powerlimits]((-1,+1))
yf[:set_scientific](true)
yf[:set_useOffset](true)
fig = figure(figsize=(8,4))
plot((0:nt)*Δt, ΔI)
ax = gca()
ax[:yaxis][:set_major_formatter](yf)
savefig("example_2nd_trapezoidal.svg")
ApproxFun
Create Poincaré Invariant
pinv = PoincareInvariant2ndApproxFun(
x₀ -> ode = ODE(f, x₀), surface, omega,
Δt, nd, n₀, n₀, nt)Compute solution
tab = TableauGLRK(1)
int = Integrator(pinv.equ, tab, pinv.Δt)
sol = integrate(pinv.equ, int, pinv.ntime)
nothing # hidePlot solution
fig = figure(figsize=(4,4))
plot(sol.q[1,0,:], sol.q[2,0,:], ".", markersize=.5)
savefig("example_2nd_approxfun_initial_state.svg")
fig = figure(figsize=(4,4))
plot(sol.q[1,end,:], sol.q[2,end,:], ".", markersize=.5)
savefig("example_2nd_approxfun_final_state.svg")
fig = figure(figsize=(4,4))
for i in 1:n₀
plot3D(sol.q[1,:,i], sol.q[2,:,i], collect(0:nt)*Δt)
end
savefig("example_2nd_approxfun_evolution1.png")
fig = figure(figsize=(4,4))
for i in 0:20:nt
plot3D(sol.q[1,i,:], sol.q[2,i,:], i*Δt)
end
savefig("example_2nd_approxfun_evolution2.png")
Compute Poincarè invariant
I, J, ΔI, ΔJ = evaluate_poincare_invariant(pinv, sol)
pinv.I[0]Plot invariant error
yf = matplotlib[:ticker][:ScalarFormatter]()
yf[:set_powerlimits]((-1,+1))
yf[:set_scientific](true)
yf[:set_useOffset](true)
fig = figure(figsize=(8,4))
plot((0:nt)*Δt, ΔI)
ax = gca()
ax[:yaxis][:set_major_formatter](yf)
savefig("example_2nd_approxfun.svg")
OrthogonalPolynomialsQuasi
Create Poincaré Invariant
pinv = PoincareInvariant2ndOPQ(
x₀ -> ode = ODE(f, x₀), surface, omega,
Δt, nd, n₀, n₀, nt)
nothing # hideCompute solution
tab = TableauGLRK(1)
int = Integrator(pinv.equ, tab, pinv.Δt)
sol = integrate(pinv.equ, int, pinv.ntime)
nothing # hidePlot solution
fig = figure(figsize=(4,4))
plot(sol.q[1,0,:], sol.q[2,0,:], ".", markersize=.5)
savefig("example_2nd_opq_initial_state.svg")
fig = figure(figsize=(4,4))
plot(sol.q[1,end,:], sol.q[2,end,:], ".", markersize=.5)
savefig("example_2nd_opq_final_state.svg")
fig = figure(figsize=(4,4))
for i in 1:n₀
plot3D(sol.q[1,:,i], sol.q[2,:,i], collect(0:nt)*Δt)
end
savefig("example_2nd_opq_evolution1.png")
fig = figure(figsize=(4,4))
for i in 0:20:nt
plot3D(sol.q[1,i,:], sol.q[2,i,:], i*Δt)
end
savefig("example_2nd_opq_evolution2.png")
Compute Poincarè invariant
I, J, ΔI, ΔJ = evaluate_poincare_invariant(pinv, sol)
pinv.I[0]Plot invariant error
yf = matplotlib[:ticker][:ScalarFormatter]()
yf[:set_powerlimits]((-1,+1))
yf[:set_scientific](true)
yf[:set_useOffset](true)
fig = figure(figsize=(8,4))
plot((0:nt)*Δt, ΔI)
ax = gca()
ax[:yaxis][:set_major_formatter](yf)
savefig("example_2nd_opq.svg")