The Discretized Linear Wave

The discretized version of the 1d linear wave equation with homogeneous Dirichlet boundary conditions.

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

The system is discretized on N equidistant interior points $ξ_i = iΔξ - 1/2$, $i = 1, …, N$, with $Δξ = 1/(N+1)$; N is the $\tilde{N}$ of the documentation page. The two boundary values $q_0 = q_{N+1} = 0$ are fixed and are not part of the state, which therefore has exactly N components. This is the full-order model of [2]: the Hamiltonian is $H = \tfrac12 p^Tp - \tfrac{μ^2}{2} q^T D_{ξξ} q$ with the central-difference matrix $D_{ξξ}$ and the canonical Poisson matrix.

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. (These figures were measured before 0.10, with the two boundary points still in the state; the scaling is unchanged.)

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₀ = initial_condition(μ̃, N).q,
    p₀ = initial_condition(μ̃, N).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.

initial_condition(μ, N)

The initial condition $q(0, ξ_i) = u_0(ξ_i; μ)$, $p(0, ξ_i) = -μ ∂_ξ u_0(ξ_i; μ)$ on the N interior_domain points, as a named tuple (q, p). Before it reaches a boundary, the solution of the continuous problem from this initial condition is the travelling wave $u_0(ξ - μt; μ)$.

interior_domain(N, T = Float64)

The N interior grid points $ξ_i = i/(N+1) - 1/2$, $i = 1, …, N$, at which the state is defined. compute_domain(N + 2) is the same grid with the two boundary points included.

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₀ = initial_condition(μ̃, N).q,
    p₀ = initial_condition(μ̃, N).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 with $q_0 = q_{N+1} = 0$,

\[V(q) = \frac{\mu^2}{2Δξ^2} \sum_{i = 0}^{N} (q_{i+1} - q_i)^2 = -\frac{\mu^2}{2} q^T D_{ξξ} q, \qquad Δξ = \frac{1}{N+1},\]

where $D_{ξξ}$ is the $N × N$ central-difference matrix $\mathrm{tridiag}(1, -2, 1)/Δξ^2$.

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.

ω!(Ω, 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$. 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 = -μ^2 D_{ξξ} q$ of the discretized gradient energy potential, written into dV. α = -1 gives the force $μ^2 D_{ξξ} q$ directly, which is what the force functions want and saves a second pass over the state compared with negating afterwards.

The boundary values $q_0 = q_{N+1} = 0$ enter the stencil at the first and the last point, so the two end rows drop one neighbor and the loop over the interior rows is branch-free. test/linear_wave.jl checks this against a central finite difference of potential at several N, including N = 1, where the two end rows coincide.