The Discretized Linear Wave

The discretized version of the 1d linear wave equation.

It is a prime example of a non-trivial completely integrable system.

The system is discretized on N interior points — the $\tilde{N}$ of the documentation page — so the state has N + 2 components once the two boundary points are included.

Like the number of lattice sites of TodaLattice, N fixes the size of the system: it sets the number of degrees of freedom and the summation bounds, so it cannot survive symbolize and hence cannot be a system parameter. It is instead a leading positional argument of the problem constructors, defaulting to Ñ = 256, and hamiltonian/lagrangian take it as a trailing argument. The only system parameter proper is $\mu$.

hodeproblem, lodeproblem, hodeensemble and lodeensemble use hand-written vector fields by default: in-place v, f, ϑ, g and ω! over a shared ∇V! kernel. Passing symbolic = true generates the equations of motion with EulerLagrange via hamiltonian_system or lagrangian_system instead, which agrees to round-off and is kept for cross-checking.

The symbolic route is unusable at the default size. EulerLagrange builds the $2n × 2n$ Lagrange two-form by symbolically differentiating a dense matrix, which at $n = 258$ is 266 k entries: lagrangian_system(256, …) takes 155 s and emits 14 MB of code for an ω whose first evaluation costs a further 147 s — all for a two-form that no integrator in GeometricIntegrators ever evaluates and that GeometricEquations.check_methods skips. Construction grows as $n^{2.5}$ and the generated ω as $n^2$, while hodeproblem/lodeproblem in their default form build in a tenth of a millisecond at any size. Per call the hand-written code is the faster one too, by 1.4–1.9× for the force and 6–8× for H and L. See benchmark/linear_wave.jl.

Produces initial conditions for the bump function. Here the $p$-part is initialized with zeros.

Produces initial condition for the bump function. Here the $p$-part is initialized as $-\mu\partial_\xi u_0(\xi, \mu)$.

Third-degree spline that is used as a basis to construct the initial conditions.

hamiltonian_system(N, parameters)

The EulerLagrange HamiltonianSystem for the linear wave equation on N interior points, from which hodeproblem(…; symbolic = true) takes its v, f and H.

Cheap at every size — 1.9 s at the default Ñ = 256 — since a Hamiltonian system carries no two-form. Contrast lagrangian_system.

hodeensemble(N, q₀, p₀; timespan, timestep, parameters, symbolic)

Hamiltonian ensemble for the linear wave equation (varying initial conditions and/or parameters).

Takes the same arguments as hodeproblem, with q₀, p₀ and/or parameters given as vectors of samples.

hodeproblem(N, q₀, p₀; timespan, timestep, parameters, symbolic)

Hamiltonian problem for the linear wave equation on N interior points.

Constructor with default arguments:

hodeproblem(
    N  = 256,
    q₀ = compute_initial_condition2(μ̃, N + 2).q,
    p₀ = compute_initial_condition2(μ̃, N + 2).p;
    timespan   = (0, 1),
    timestep   = 0.005025125628140704,
    parameters = (μ = 0.6,),
    symbolic   = false
)

With symbolic = true the equations of motion are generated with EulerLagrange via hamiltonian_system instead of using the hand-written vector fields. The two agree to round-off; the symbolic route is kept for cross-checking and costs 1.9 s to build at N = 256.

lagrangian_system(N, parameters)

The EulerLagrange LagrangianSystem for the linear wave equation on N interior points, from which lodeproblem(…; symbolic = true) takes its ϑ, f, g, ω and L.

This is the expensive one: it builds the $2n × 2n$ two-form by differentiating a dense matrix, so its cost grows as $n^{2.5}$ and reaches 155 s at the default Ñ = 256. See the module docstring and benchmark/linear_wave.jl.

lodeensemble(N, q₀, p₀; timespan, timestep, parameters, symbolic)

Lagrangian ensemble for the linear wave equation (varying initial conditions and/or parameters).

Takes the same arguments as lodeproblem, with q₀, p₀ and/or parameters given as vectors of samples.

lodeproblem(N, q₀, p₀; timespan, timestep, parameters, symbolic)

Lagrangian problem for the linear wave equation on N interior points.

Constructor with default arguments:

lodeproblem(
    N  = 256,
    q₀ = compute_initial_condition2(μ̃, N + 2).q,
    p₀ = compute_initial_condition2(μ̃, N + 2).p;
    timespan   = (0, 1),
    timestep   = 0.005025125628140704,
    parameters = (μ = 0.6,),
    symbolic   = false
)

With symbolic = true the equations of motion are generated with EulerLagrange via lagrangian_system instead of using the hand-written vector fields. The two agree to round-off, but the symbolic route takes 155 s to build at N = 256 — see the module docstring.

potential(q, parameters, N)

The discretized gradient energy on N interior points,

\[V(q) = \frac{\mu^2}{4Δξ^2} \sum_{i = 2}^{N+1} \left[ (q_i - q_{i-1})^2 + (q_{i+1} - q_i)^2 \right], \qquad Δξ = \frac{1}{N+1}.\]

Both the Hamiltonian and the Lagrangian take N as a trailing argument (as in TodaLattice) rather than reading it from parameters: it fixes the number of degrees of freedom and the summation bounds, so it has to be a plain integer and cannot survive symbolize.

The sum runs over the interior points only, which is what gives the two boundary points a different stencil weight from the interior — see ∇V!.

ω!(Ω, t, q, params)

The Lagrange two-form $ω = dθ$ of this regular Lagrangian, where $θ = ϑ_i \, dq^i$.

A regular Lagrangian is a second-order system of $n$ equations, equivalent to a first-order system of $2n$, so Ω is the $2n × 2n$ form on $(q, \dot q)$ with $n = N + 2$. In block form it is $[∂ϑ_i/∂q_j - ∂ϑ_j/∂q_i \; -M^T; \; M \; 0]$; since $ϑ = ∂L/∂\dot q = \dot q$ depends on the velocities alone and the mass matrix is the identity, the upper-left block vanishes and what remains is the canonical $[0 \; -I; \; I \; 0]$.

This is the convention EulerLagrange's LagrangianSystem produces, so lodeproblem agrees whether or not symbolic = true. Note that no integrator in GeometricIntegrators evaluates ω and GeometricEquations.check_methods skips it — which is what makes the 14 MB of code EulerLagrange emits for it at the default size, and the 147 s its first evaluation then costs, pure overhead.

∇V!(dV, q, parameters, N, α = 1)

$α$ times the gradient $∂V/∂q$ of the discretized gradient energy potential, written into dV. α = -1 gives the force $-∂V/∂q$ directly, which is what the force functions want and saves a second pass over the state compared with negating afterwards.

Writing $d_k = q_{k+1} - q_k$ and $c = μ^2 / 2Δx^2$, the sum in $V$ runs over the interior points $i = 2, …, N+1$ and therefore counts every interior difference twice, so $V = \tfrac{c}{2} \sum_{k=1}^{N+1} w_k d_k^2$ with $w_k = 2 - [k = 1] - [k = N+1]$: the two boundary weights are 1, not 2. Hence $∂V/∂q_j = c \, (w_{j-1} d_{j-1} - w_j d_j)$, with $d_0 = d_{N+2} = 0$.

Substituting the weights splits that into a uniform weight-two stencil plus four scalar corrections at the boundaries. Written that way the loop is branch-free and vectorizes, each output is written exactly once, and — unlike a peeled-boundary stencil — the expression is correct for every $N ≥ 1$ with no case analysis. It is also about twice as fast as accumulating difference by difference (dV[i] += a; dV[i-1] -= a), which reads and writes adjacent elements and so cannot be vectorized, and about twice as fast as the code EulerLagrange generates.

The four correction signs are easy to get wrong and are not checked by anything else, so test/linear_wave_tests.jl verifies this against a central finite difference of potential at several N, including the N = 1 case that a peeled-boundary stencil gets wrong.