Charged Particle

In the second part of the tutorial we will treat a non-canonical system and show how to compute! Poincaré integral invariants for it. The system consists of a 2D charged particle subject to electromagnetic fields. The equations of motion are

\[\dot{x} = v,\; \dot{v} = E(x) + v × B(x)\]

where $x$ is the position vector, $v$ is the velocity, $E$ and $B$ are the electric and magnetic field and the mass and charge of the particle is set to one. For our example we will use a constant magnetic field of strength $10$ pointing in the z-direction and an electric field given by $E(x, y) = (-x, -y^3)$.

Integration

Again, we'll quickly implement our own integrators here, namely the Runge-Kutta method (RK4) and a second order symplectic splitting method for charged particles in electromagnetic fields.

using PoincareInvariants
using StaticArrays

E(x, y) = (-x, -y^3)

A(x, y) = 5 .* (-y, x)
B(x, y) = 10.0
Bx(x, y, Δx) = 10 * Δx  # integral in x direction from x to x + Δx
By(x, y, Δy) = 10 * Δy  # integral in y direction from y to y + Δy

struct Split2 end

function ϕx((x, y, vx, vy), dt)
    Δx = vx * dt
    return (x + Δx, y, vx, vy - Bx(x, y, Δx))
end

function ϕy((x, y, vx, vy), dt)
    Δy = vy * dt
    return (x, y + Δy, vx + By(x, y, Δy), vy)
end

function ϕE((x, y, vx, vy), dt)
    (Δvx, Δvy) = dt .* E(x, y)
    return (x, y, vx + Δvx, vy + Δvy)
end

function timestep(z, dt, ::Split2)
    hdt = 0.5 * dt
    z = ϕx(z, hdt)
    z = ϕy(z, hdt)
    z = ϕE(z, dt)
    z = ϕy(z, hdt)
    z = ϕx(z, hdt)
    return z
end

struct RK4 end

function zdot((x, y, vx, vy))
    ex, ey = E(x, y); b = B(x, y)
    (vx, vy, ex + b * vy, ey - b * vx)
end

function timestep(z, dt, ::RK4)
    hdt = 0.5 .* dt
    k1 = zdot(z)
    k2 = zdot(z .+ hdt .* k1)
    k3 = zdot(z .+ hdt .* k2)
    k4 = zdot(z .+  dt .* k3)
    return z .+ dt .* (k1 .+ 2 .* k2 .+ 2 .* k3 .+ k4) .* (1/6)
end

"""
    integrate(z0, dt, nsteps, nt, method)

start at `z0` and integrate the equations of motion using `method`.
Returns the timeseries as a vector of tuples. `nt` points are saved,
`nsteps` steps are taken from saved point to saved point and `dt` is the
size of each time step.
"""
function integrate(z0, dt, nsteps, nt, method)
    out = Vector{NTuple{4, Float64}}(undef, nt)
    z = out[1] = z0
    for i in 2:nt
        for _ in 1:nsteps
            z = timestep(z, dt, method)
        end
        out[i] = z
    end
    return out
end

integrate(mat::AbstractMatrix, dt, nsteps, nt, method) = map(eachrow(mat)) do r
    integrate((r[1], r[2], r[3], r[4]), dt, nsteps, nt, method)
end

Computing the Invariants

Having gotten integration out of the way, we can move onto the invariants. For this system, the first invariant is given by

\[I_{1} = \int_{\gamma} v(x) + A(x) \, dx\]

where $A(x)$ is the magnetic vector potential as a function of the position vector $x$. The second invariant is given by

\[I_{2} = \int_{S} \omega_{ij} (q) \, dz^{i} \, dz^{j}\]

with the two-form $\omega$ given by

\[\begin{pmatrix} 0 & B & -1 & 0 \\ -B & 0 & 0 & -1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{pmatrix}\]

In code, we can define these forms as follows.

function oneform((x, y, vx, vy), ::Real, ::Any)
    p = (vx, vy) .+ A(x, y)
    @SVector [p[1], p[2], 0, 0]
end

function twoform(z, ::Real, ::Any)
    b = B(z[1], z[2])
    @SMatrix [ 0  b -1  0;
              -b  0  0 -1;
               1  0  0  0;
               0  1  0  0]
end

Forms always have the signature form(z, t, p), where z is the phasespace position, t is a time and p is an arbitrary parameter. To use these forms to compute our invariants, we must create the setup objects FirstPI and SecondPI and initialise the curve and surface we want in phasespace.

pi1 = FirstPI{Float64, 4}(oneform, 1_000)
pi2 = SecondPI{Float64, 4}(twoform, 10_000)

I1 = 0
pnts1 = getpoints(pi1) do θ
    10 .* (0, 0, sinpi(2θ), cospi(2θ))
end

I2 = 1_000
pnts2 = getpoints(pi2) do x, y
    10 .* (y - 0.5, x - 0.5, 0, 0)
end

The curve is just a point in space, so the first integral invariant evaluates to zero. The second invariant is equal to $1000$ since, it is equal to the magnetic field of strength $10$ integrated over a $10$ by $10$ square. We can confirm this by computing the inital invariants.

@assert isapprox(compute!(pi1, pnts1, 52.3, "optional parameter"), I1; atol=10^(-15))
@assert isapprox(compute!(pi2, pnts2), I2; atol=10^(-11))

By default p is nothing and t is zero, but as shown in the first line we can also explicitly supply a time and parameter if the forms we defined earlier had required them. compute! used on a time series accepts an iterable of times and an arbitrary parameter.

times = range(0.0; step=0.05 * 50, length=5)

series1 = integrate(pnts1, 0.05, 50, 5, Split2())
@assert all(compute!(pi1, series1, times, ("optional parameters", 3.7, 42))) do I
    abs(I - I1) < 10^(-13)
end

series2 = integrate(pnts2, 0.05, 50, 5, Split2())
@assert all(compute!(pi2, series2)) do I
    abs(I - I2) < 10^(-11)
end

Again, we can plot the results. The top two plots show the increasingly distorted curve and surface over time (projected onto x and y position components), while the bottom two plots show the error in the invariant over time for the two integration algorithms.

We see that RK4 does not preserve the non-canonical invariant while the splitting method does.