Legendre-Galerkin Discretization

• Legendre-Galerkin discretization is a spectral Galerkin method that uses finite-dimensional polynomial spaces based on Legendre polynomials to enforce boundary or initial conditions exactly.
• It produces structured operator matrices—often sparse—with engineered derivative properties that enable efficient, adaptive, and conservation-preserving numerical solutions.
• The methodology is applied in variational settings for elliptic and time-dependent PDEs, offering robust frameworks for handling nonlinearities, preconditioning, and adaptive refinement.

Searching arXiv for recent and foundational papers on Legendre-Galerkin discretization to ground the article. Legendre-Galerkin discretization is a spectral Galerkin methodology in which the trial and test spaces are finite-dimensional polynomial spaces generated by Legendre polynomials or by modified Legendre combinations chosen to satisfy boundary or initial conditions exactly. In the formulations represented across current work, the continuous problem is first written in weak form and then restricted to spaces such as $V_N\subset P_N$, tensor-product polynomial spaces on $(-1,1)^d$, or Petrov-Galerkin pairs tailored to initial value problems; the resulting semidiscrete or fully discrete systems are governed by mass, stiffness, and differentiation matrices whose structure depends strongly on the basis design (Kong et al., 2022, Diao et al., 2020, Liu et al., 2016). A central theme is that Legendre-based bases can be engineered so that derivative operators become diagonal or even the identity, while at the same time preserving exact boundary constraints, enabling sparse or structured linear systems, and supporting adaptive, conservation-preserving, and data-driven variants (Kong et al., 2022, Canuto et al., 2012, Ko et al., 2022).

1. Variational setting and finite-dimensional polynomial spaces

The common starting point is a weak formulation posed on a Hilbert space adapted to the PDE and its constraints. For second-order elliptic problems on $I=(-1,1)$, one representative form is

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$

with $H=H_0^1(I)$ in the Dirichlet case and a suitable $H^1(I)$-based space in the Neumann case (Ko et al., 2022). For self-adjoint elliptic problems on $\Omega=(-1,1)^d$,

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$

and the Legendre-Galerkin approximation seeks $u_N\in X_N^d$ satisfying the same identity for all $v\in X_N^d$ (Diao et al., 2020). For time-dependent problems, the same principle appears either in semidiscrete form,

$(-1,1)^d$0

or as a space-time formulation in which time itself is discretized spectrally (Płociniczak, 2021, Kong et al., 2022).

A general Galerkin truncation may be written as

$(-1,1)^d$1

with mass matrix $(-1,1)^d$2 and reduced system

$(-1,1)^d$3

This coefficient-space representation is the basic algebraic form into which Legendre-Galerkin discretizations fit naturally when $(-1,1)^d$4 are Legendre polynomials or modified Legendre functions satisfying the required constraints (Hilliard et al., 2023).

In initial value problems, the Legendre dual-Petrov-Galerkin construction uses different trial and test spaces. For the first-order model problem

$(-1,1)^d$5

the approximation is written as

$(-1,1)^d$6

with dual test space

$(-1,1)^d$7

and Petrov-Galerkin formulation

$(-1,1)^d$8

This distinguishes the Legendre dual-Petrov-Galerkin framework from conforming symmetric Galerkin formulations used for elliptic PDEs (Kong et al., 2022).

2. Basis construction and exact enforcement of constraints

A defining feature of Legendre-Galerkin discretization is the use of modified Legendre bases that encode the boundary or initial conditions directly. For homogeneous Dirichlet conditions on $(-1,1)^d$9, a standard choice is

$I=(-1,1)$0

which satisfies $I=(-1,1)$1 because $I=(-1,1)$2; this yields

$I=(-1,1)$3

and extends by tensor products in two and three dimensions (Diao et al., 2020). A normalized Shen-type variant used for the two-dimensional Schrödinger equation is

$I=(-1,1)$4

again chosen so that $I=(-1,1)$5 (Liu et al., 2016).

Other problems motivate different Legendre subspaces. For the degenerate climatology model, the approximation space is

$I=(-1,1)$6

with basis

$I=(-1,1)$7

so that $I=(-1,1)$8 and

$I=(-1,1)$9

The convenience of this choice follows directly from the degenerately weighted operator (Płociniczak, 2021).

For 1D elliptic Dirichlet and Neumann problems, another widely used modified basis is

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$0

with coefficients selected to satisfy the boundary conditions exactly. The explicit choices are

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$1

for the Dirichlet case, and

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$2

for the Neumann case (Ko et al., 2022).

In adaptive Legendre-Galerkin theory, the Babuška-Shen basis plays an analogous role for $$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$3: $$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$4 with

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$5

This gives an $$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$6-orthonormal modal basis in one dimension and underlies coefficient-space residual estimates and adaptive algorithms (Canuto et al., 2012).

The initial-value setting admits yet another specialized basis design. For the first-order dual-Petrov-Galerkin scheme,

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$7

which produces the relation

$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$8

This basis choice is not merely a convenience; it is the mechanism by which the highest derivative matrix becomes the identity (Kong et al., 2022).

3. Matrix structures, sparsity, and spectral algebra

The algebraic structure of a Legendre-Galerkin discretization is determined by the interaction between the PDE operator and the selected basis. Several recurrent patterns appear.

Setting	Legendre-based basis	Induced matrix structure
First-order LDPG IVP	$$\int_I \left(u'(x)v'(x)+\nu u(x)v(x)\right)\,dx = \int_I f(x)v(x)\,dx \qquad \forall v\in H,$$9 with $H=H_0^1(I)$0	derivative matrix $H=H_0^1(I)$1, mass matrix tridiagonal
1D elliptic Dirichlet/Neumann	$H=H_0^1(I)$2	stiffness matrix diagonal, mass matrix symmetric penta-diagonal
2D Schrödinger spatial discretization	$H=H_0^1(I)$3 tensor products	derivative inner products diagonal, tensor-product mass/stiffness form
Variable-coefficient elliptic problem	$H=H_0^1(I)$4 tensor products	exact Galerkin matrix dense and ill-conditioned
Adaptive $H=H_0^1(I)$5 setting	Babuška-Shen basis $H=H_0^1(I)$6	coefficient norm representation; sparse mass matrix in 1D

In the first-order LDPG scheme, writing

$H=H_0^1(I)$7

yields

$H=H_0^1(I)$8

The matrix $H=H_0^1(I)$9 is tridiagonal and nonsymmetric, and $H^1(I)$0 is exactly the Jacobi matrix associated with the three-term recurrence of the generalized Bessel polynomials $H^1(I)$1 for $H^1(I)$2 and $H^1(I)$3. Consequently, the eigenvalues of $H^1(I)$4 are exactly the negatives of the zeros of $H^1(I)$5, the eigenvectors are obtained by polynomial evaluation, and the eigenvalues are simple, conjugate-symmetric, and located in the open right half-plane (Kong et al., 2022). For second-order IVPs, the analogous matrix is penta-diagonal and satisfies

$H^1(I)$6

where $H^1(I)$7 is the Jacobi matrix for $H^1(I)$8 (Kong et al., 2022).

For elliptic problems discretized with Legendre combinations satisfying boundary constraints, sparsity can be equally pronounced. In the ULGNet discretization of $H^1(I)$9, the discrete system is

$\Omega=(-1,1)^d$0

with

$\Omega=(-1,1)^d$1

Here the stiffness matrix is diagonal,

$\Omega=(-1,1)^d$2

while the mass matrix is symmetric penta-diagonal with explicitly listed nonzero bands (Ko et al., 2022). In the two-dimensional Schrödinger problem, the normalized Shen basis satisfies

$\Omega=(-1,1)^d$3

while the one-dimensional mass matrix $\Omega=(-1,1)^d$4 is sparse and tridiagonal-like (Liu et al., 2016).

These sparse structures are not universal. For the non-separable elliptic equation

$\Omega=(-1,1)^d$5

the exact Legendre-Galerkin matrix $\Omega=(-1,1)^d$6 is dense and ill-conditioned when the coefficients are variable and non-separable (Diao et al., 2020). This is a recurring correction to the common assumption that Legendre orthogonality alone guarantees simple linear algebra. A plausible implication is that basis orthogonality controls only part of the operator representation; variable coefficients can still destroy sparsity unless additional approximation or preconditioning is introduced.

4. Time-dependent formulations and space-time discretization

Legendre-Galerkin discretization is used in both semidiscrete-in-space and fully space-time settings. In the two-dimensional Schrödinger equation, after mapping the physical square to $\Omega=(-1,1)^d$7 and subtracting a boundary-lifting function, the homogeneous problem is discretized in the tensor-product space

$\Omega=(-1,1)^d$8

The expansion

$\Omega=(-1,1)^d$9

produces the matrix ODE

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$0

or, after vectorization,

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$1

Time integration is then performed by an implicit Runge-Kutta method (Liu et al., 2016).

The LDPG framework goes further by discretizing time spectrally. For a linear evolutionary system

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$2

the time approximation

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$3

leads to

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$4

Two solvers are emphasized. Matrix diagonalization uses $$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$5 and decouples the problem into

$$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$6

This is fully parallel, but $$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$7 is extremely ill-conditioned and is practical only for modest $$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$8. The alternative is a QZ or generalized Schur decomposition, which is more stable for large $$\int_\Omega B(x)\nabla u\cdot\nabla v\,dx+\int_\Omega a(x)uv\,dx=\int_\Omega fv\,dx \qquad \forall v\in H_0^1(\Omega),$$9 because it avoids explicit inversion and reduces the system by backward substitution (Kong et al., 2022).

For nonlinear and nonlocal parabolic equations, the spatial Legendre-Galerkin discretization is combined with a weighted $u_N\in X_N^d$0-scheme, extrapolation, and quadrature for the memory term. The fully discrete method is

$u_N\in X_N^d$1

with

$u_N\in X_N^d$2

Here $u_N\in X_N^d$3 gives a Crank-Nicolson-type scheme with second-order time accuracy, and the extrapolation makes the scheme linear at each time step (Płociniczak, 2021).

5. Elliptic operators, preconditioning, and adaptive Legendre-Galerkin methods

For non-separable variable-coefficient elliptic operators, the principal difficulty is not the polynomial approximation itself but the linear system generated by it. The preconditioned Legendre spectral Galerkin method addresses

$u_N\in X_N^d$4

on $u_N\in X_N^d$5, $u_N\in X_N^d$6, with the Dirichlet-adapted basis $u_N\in X_N^d$7. The exact system is symmetric positive definite but dense and ill-conditioned. The proposed preconditioner $u_N\in X_N^d$8 is built by approximating $u_N\in X_N^d$9 and $v\in X_N^d$0 with truncated Legendre series, producing a banded or block-banded matrix whose bandwidth depends on the truncation degrees rather than on $v\in X_N^d$1. Fast matrix-vector multiplication is performed without explicit formation of the dense matrix by using backward discrete Legendre transforms, pointwise multiplication at Legendre-Gauss nodes, forward discrete Legendre transforms, and short recurrences. The resulting complexity is

$v\in X_N^d$2

for operator application, while the ILU(0)-based approximate preconditioner solve costs

$v\in X_N^d$3

and the total PCG complexity is also

$v\in X_N^d$4

for $v\in X_N^d$5 (Diao et al., 2020).

A distinct development concerns adaptive Legendre-Galerkin methods. In one dimension, the Babuška-Shen basis gives a coefficient representation of the $v\in X_N^d$6-norm and permits residual-based adaptivity for

$v\in X_N^d$7

The Galerkin approximation on

$v\in X_N^d$8

satisfies residual-error equivalence, and the ideal adaptive algorithm based on Dörfler marking yields the contraction

$v\in X_N^d$9

A predictor-corrector version augments the marked set by enrichment and then applies coarsening; it also satisfies linear convergence and is proved optimal in cardinality relative to Gevrey-type sparsity classes (Canuto et al., 2012).

In multiple dimensions, tensorized Babuška-Shen functions are no longer directly usable as an $(-1,1)^d$00-orthogonal basis. The multidimensional adaptive theory therefore introduces a quasi-orthonormalization that constructs a nearly-orthonormal Babuška-Shen basis $(-1,1)^d$01 from the tensorized basis. This yields norm equivalences

$(-1,1)^d$02

which enable feasible residual computation, Dörfler marking, enrichment, coarsening, and contraction: $(-1,1)^d$03 This places multidimensional adaptive Legendre-Galerkin discretization within the framework of nonlinear approximation and Gevrey-type sparsity (Canuto et al., 2014).

6. Conservation, stability pathologies, and related hybrid formulations

One persistent issue is that standard Galerkin truncation does not in general preserve the first integrals of the underlying PDE. In the general reduced system

$(-1,1)^d$04

the Galerkin RONS modification enforces invariants $(-1,1)^d$05 by solving the constrained problem

$(-1,1)^d$06

which yields

$(-1,1)^d$07

The correction term is chosen so that

$(-1,1)^d$08

holds at the continuous-time reduced level. Since the paper does not specialize to Legendre polynomials, this is a general Galerkin framework, but it applies directly when the basis $(-1,1)^d$09 is Legendre or Legendre-like (Hilliard et al., 2023).

Stability questions can also be more subtle than approximation estimates alone suggest. For the 1D homogeneous wave equation discretized by a classical Legendre-Galerkin semidiscretization, the high-frequency numerical modes are spurious from the point of view of wave propagation and destroy uniform boundary observability. The paper states that, for the standard scheme,

$(-1,1)^d$10

Three remedies are studied numerically: spectral filtering, a mixed formulation, and Nitsche’s method. Each recovers a uniform positive lower bound on the observability constant, but none appears to restore a uniform direct or trace inequality (Gagnon et al., 2016). This directly counters the misconception that spectral accuracy alone guarantees stable control-theoretic behavior.

A further line of development uses Legendre-Galerkin discretization as the finite-dimensional backbone of neural methods. In unsupervised Legendre-Galerkin neural networks, the solution is expanded as

$(-1,1)^d$11

and the network is trained to predict $(-1,1)^d$12 by minimizing the residual of the discrete linear system

$(-1,1)^d$13

The population loss is exactly the squared algebraic residual, and the analysis proves that the empirical minimizer converges to the Legendre-Galerkin weak solution $(-1,1)^d$14 under the stated assumptions (Ko et al., 2022). This suggests that, in this setting, the neural component parameterizes the coefficient map rather than replacing the underlying spectral discretization.

A related but distinct methodology appears in discontinuous Galerkin discretizations based on hierarchical Legendre polynomial basis functions. In the $(-1,1)^d$15-version SIPDG setting, the Legendre basis is modal, orthogonal in $(-1,1)^d$16, and hierarchical, so prolongation between polynomial levels is the natural injection and the mass matrix is diagonal. The paper explicitly characterizes this as “Legendre-Galerkin-like” while emphasizing that the formulation remains discontinuous and element-local rather than conforming (Lei et al., 17 Sep 2025). This establishes a clear boundary of the topic: hierarchical Legendre modal structure can persist beyond conforming Galerkin methods, but the continuity constraints and bilinear forms change fundamentally.

Legendre-Galerkin discretization therefore encompasses more than a single scheme. Across current formulations it denotes a family of modal polynomial Galerkin constructions in which Legendre-based trial and test spaces are adapted to operator structure, constraints, and computational objectives. Depending on the problem, the method may deliver exact derivative identities, sparse or structured operator matrices, spectral convergence, adaptive contraction, conservation corrections, or stable space-time solvers; it may also exhibit dense matrices, ill-conditioning, or high-frequency pathologies that require preconditioning, filtering, or reformulation (Kong et al., 2022, Diao et al., 2020, Gagnon et al., 2016).