---
title: Legendre-Galerkin Discretization
url: https://www.emergentmind.com/topics/legendre-galerkin-discretization
type: topic
---

# Legendre-Galerkin Discretization

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 [2211.10599], [2004.13961], [1607.08420]. 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 [2211.10599], [1206.5524], [2211.08900].

## 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 [2211.08900]. 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\) [2004.13961]. For time-dependent problems, the same principle appears either in semidiscrete form,
\[
(u_{N,t},v)+a(D(u_N);u_N,v)=(f(t,u_N,Ju_N),v)\qquad \forall v\in V_N,
\]
or as a space-time formulation in which time itself is discretized spectrally [2106.05140], [2211.10599].

A general Galerkin truncation may be written as
\[
\hat u(x,a(t))=\sum_{i=1}^N a_i(t)\phi_i(x),
\]
with mass matrix \(M_{ij}=\langle \phi_i,\phi_j\rangle_H\) and reduced system
\[
M\dot a=f(a).
\]
This coefficient-space representation is the basic algebraic form into which Legendre-Galerkin discretizations fit naturally when \(\phi_i\) are Legendre polynomials or modified Legendre functions satisfying the required constraints [2312.05322].

In initial value problems, the Legendre dual-Petrov-Galerkin construction uses different trial and test spaces. For the first-order model problem
\[
u'(t)=\sigma u(t),\qquad u(-1)=u_0,
\]
the approximation is written as
\[
u_N(t)=u_0+v_N(t),\qquad v_N\in {}_0\mathbb P_N:=\{\phi\in\mathbb P_N:\phi(-1)=0\},
\]
with dual test space
\[
{}^0\mathbb P_N:=\{\psi\in\mathbb P_N:\psi(1)=0\},
\]
and Petrov-Galerkin formulation
\[
(v_N',\psi)-\sigma(v_N,\psi)=\sigma(u_0,\psi),\qquad \forall \psi\in {}^0\mathbb P_N.
\]
This distinguishes the Legendre dual-Petrov-Galerkin framework from conforming symmetric Galerkin formulations used for elliptic PDEs [2211.10599].

## 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)\), a standard choice is
\[
\phi_k(x)=L_k(x)-L_{k+2}(x),\qquad k=0,1,\dots,N-2,
\]
which satisfies \(\phi_k(\pm1)=0\) because \(L_n(\pm 1)=(\pm1)^n\); this yields
\[
X_N=\{v\in P_N:\ v(\pm1)=0\}=\operatorname{span}\{\phi_0,\dots,\phi_{N-2}\}
\]
and extends by tensor products in two and three dimensions [2004.13961]. A normalized Shen-type variant used for the two-dimensional Schrödinger equation is
\[
\phi_k(\tilde x)=c_k\big(L_k(\tilde x)-L_{k+2}(\tilde x)\big), \qquad c_k=\frac{1}{\sqrt{4k+6}},
\]
again chosen so that \(\phi_k(\pm1)=0\) [1607.08420].

Other problems motivate different Legendre subspaces. For the degenerate climatology model, the approximation space is
\[
V_N=\{v\in \mathbb{P}_N(0,1): v_x(0)=0,\ v_x(1)<\infty\},
\]
with basis
\[
\phi_i=\sqrt{4i+1}\,L_{2i}(x),\qquad i=0,1,2,\dots,N,
\]
so that \((\phi_i,\phi_j)=\delta_{ij}\) and
\[
L\phi_i:=-\big((1-x^2)\phi_{i,x}\big)_x=\lambda_i\phi_i,\qquad \lambda_i=2i(2i+1).
\]
The convenience of this choice follows directly from the degenerately weighted operator [2106.05140].

For 1D elliptic Dirichlet and Neumann problems, another widely used modified basis is
\[
\phi_k(x)=L_k(x)+a_k L_{k+1}(x)+b_k L_{k+2}(x),
\]
with coefficients selected to satisfy the boundary conditions exactly. The explicit choices are
\[
a_k=0,\qquad b_k=-1
\]
for the Dirichlet case, and
\[
a_k=0,\qquad b_k=\frac{-k(k+1)}{(k+2)(k+3)}
\]
for the Neumann case [2211.08900].

In adaptive Legendre-Galerkin theory, the Babuška-Shen basis plays an analogous role for \(H_0^1(I)\):
\[
\eta_k(x)=\sqrt{k-\tfrac12}\int_x^1 L_{k-1}(s)\,ds
=\frac{1}{\sqrt{4k-2}}\big(L_{k-2}(x)-L_k(x)\big),\qquad k\ge 2,
\]
with
\[
D\eta_k=-\phi_{k-1},
\qquad
(\eta_k,\eta_m)_{H_0^1(I)}=\delta_{km}.
\]
This gives an \(H_0^1\)-orthonormal modal basis in one dimension and underlies coefficient-space residual estimates and adaptive algorithms [1206.5524].

The initial-value setting admits yet another specialized basis design. For the first-order dual-Petrov-Galerkin scheme,
\[
\phi_k(t)=\frac{k+1}{\sqrt2}\big(P_k(t)+P_{k+1}(t)\big),\qquad
\phi_j^*(t)=\frac{1}{\sqrt2(j+1)}\big(P_j(t)-P_{j+1}(t)\big),
\]
which produces the relation
\[
(\phi_k',\phi_j^*)=\delta_{jk}.
\]
This basis choice is not merely a convenience; it is the mechanism by which the highest derivative matrix becomes the identity [2211.10599].

## 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 | \(\phi_k,\phi_j^*\) with \((\phi_k',\phi_j^*)=\delta_{jk}\) | derivative matrix \(=I\), mass matrix tridiagonal |
| 1D elliptic Dirichlet/Neumann | \(\phi_k=L_k+a_kL_{k+1}+b_kL_{k+2}\) | stiffness matrix diagonal, mass matrix symmetric penta-diagonal |
| 2D Schrödinger spatial discretization | \(c_k(L_k-L_{k+2})\) tensor products | derivative inner products diagonal, tensor-product mass/stiffness form |
| Variable-coefficient elliptic problem | \(\phi_k=L_k-L_{k+2}\) tensor products | exact Galerkin matrix dense and ill-conditioned |
| Adaptive \(H_0^1\) setting | Babuška-Shen basis \(\eta_k\) | coefficient norm representation; sparse mass matrix in 1D |

In the first-order LDPG scheme, writing
\[
v_N(t)=\sum_{k=0}^{N-1}\tilde v_k\phi_k(t)
\]
yields
\[
(I_N-\sigma M)\tilde{\mathbf v}=\sqrt2\,\sigma u_0\,e_1,
\qquad
M_{jk}=(\phi_k,\phi_j^*).
\]
The matrix \(M\) is tridiagonal and nonsymmetric, and \(-M\) is exactly the Jacobi matrix associated with the three-term recurrence of the generalized Bessel polynomials \(B_n^{(\alpha)}\) for \(\alpha=3\) and \(\beta=2\). Consequently, the eigenvalues of \(M\) are exactly the negatives of the zeros of \(B_N^{(3)}\), the eigenvectors are obtained by polynomial evaluation, and the eigenvalues are simple, conjugate-symmetric, and located in the open right half-plane [2211.10599]. For second-order IVPs, the analogous matrix is penta-diagonal and satisfies
\[
M^{(2)}=\widetilde M\,\widetilde M,
\]
where \(\widetilde M\) is the Jacobi matrix for \(B_n^{(4)}\) [2211.10599].

For elliptic problems discretized with Legendre combinations satisfying boundary constraints, sparsity can be equally pronounced. In the ULGNet discretization of \(-u''+\nu u=f\), the discrete system is
\[
(S+\nu M)\alpha=F,
\]
with
\[
S_{ij}=\int_I \phi_i'(x)\phi_j'(x)\,dx,\qquad
M_{ij}=\int_I \phi_i(x)\phi_j(x)\,dx.
\]
Here the stiffness matrix is diagonal,
\[
S_{kk}=-(4k+6)b_k,
\]
while the mass matrix is symmetric penta-diagonal with explicitly listed nonzero bands [2211.08900]. In the two-dimensional Schrödinger problem, the normalized Shen basis satisfies
\[
\tilde a_{jk}=\int_{-1}^1 \phi_k'(\tilde x)\phi_j'(\tilde x)\,d\tilde x = \delta_{jk},
\]
while the one-dimensional mass matrix \(\widetilde B\) is sparse and tridiagonal-like [1607.08420].

These sparse structures are not universal. For the non-separable elliptic equation
\[
-\nabla\cdot\big(B(x)\nabla u\big)+a(x)u=f,\qquad u|_{\partial\Omega}=0,
\]
the exact Legendre-Galerkin matrix \(A+B\) is dense and ill-conditioned when the coefficients are variable and non-separable [2004.13961]. 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 \((-1,1)^2\) and subtracting a boundary-lifting function, the homogeneous problem is discretized in the tensor-product space
\[
V_N^2=\operatorname{span}\{\phi_k(\tilde x)\phi_j(\tilde y):\ 0\le k,j\le N-2\}.
\]
The expansion
\[
\hat u_N(\tilde x,\tilde y,t)=\sum_{k,j=0}^{N-2}\alpha_{kj}(t)\,\phi_k(\tilde x)\phi_j(\tilde y)
\]
produces the matrix ODE
\[
-\mathrm{i}\,\widetilde B\,\alpha'(t)\,\widetilde B
+\gamma\big(\widetilde B\,\alpha(t)\,I_{N-1}+I_{N-1}\,\alpha(t)\,\widetilde B\big)-W_\alpha-\widehat F(t)=0,
\]
or, after vectorization,
\[
-\mathrm{i}\,(\widetilde B\otimes \widetilde B)\,\frac{d}{dt}\mathrm{vec}(\alpha)
+\Big(\gamma\,\widetilde B\otimes I+\gamma\,I\otimes \widetilde B-\widetilde W\Big)\mathrm{vec}(\alpha)
= \mathrm{vec}(\widehat F(t)).
\]
Time integration is then performed by an implicit Runge-Kutta method [1607.08420].

The LDPG framework goes further by discretizing time spectrally. For a linear evolutionary system
\[
\mathbf u'(t)+A\mathbf u(t)=\mathbf f(t),\qquad \mathbf u(-1)=\mathbf u_0,
\]
the time approximation
\[
\mathbf u(t)\approx \mathbf u_0+\sum_{k=0}^{N_t-1}\hat{\mathbf u}_k\phi_k(t)
\]
leads to
\[
U + M_t\,U\,A^\top = F.
\]
Two solvers are emphasized. Matrix diagonalization uses \(M_t=E\Lambda E^{-1}\) and decouples the problem into
\[
(I-\lambda_j A)\,w_j=g_j.
\]
This is fully parallel, but \(E\) is extremely ill-conditioned and is practical only for modest \(N_t\). The alternative is a QZ or generalized Schur decomposition, which is more stable for large \(N_t\) because it avoids explicit inversion and reduces the system by backward substitution [2211.10599].

For nonlinear and nonlocal parabolic equations, the spatial Legendre-Galerkin discretization is combined with a weighted \(\theta\)-scheme, extrapolation, and quadrature for the memory term. The fully discrete method is
\[
(\delta U^n,v)+a(D(\overline U^{n-\theta});\widehat U^{n-\theta},v)
=(f_h(\overline U^{n-\theta}),v),\qquad \forall v\in V_N,
\]
with
\[
\widehat U^{n-\theta}=\theta U^{n-1}+(1-\theta)U^n,\qquad
\overline U^{n-\theta}=(2-\theta)U^{n-1}-(1-\theta)U^{n-2}.
\]
Here \(\theta=\frac12\) gives a Crank-Nicolson-type scheme with second-order time accuracy, and the extrapolation makes the scheme linear at each time step [2106.05140].

## 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
\[
-\nabla\cdot\big(B(x)\nabla u\big)+a(x)u=f
\]
on \((-1,1)^d\), \(d=1,2,3\), with the Dirichlet-adapted basis \(\phi_k=L_k-L_{k+2}\). The exact system is symmetric positive definite but dense and ill-conditioned. The proposed preconditioner \(M\) is built by approximating \(B(x)\) and \(a(x)\) with truncated Legendre series, producing a banded or block-banded matrix whose bandwidth depends on the truncation degrees rather than on \(N\). 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
\[
O\!\big(N^d(\log N)^2\big)
\]
for operator application, while the ILU(0)-based approximate preconditioner solve costs
\[
O(T^{2d}N^d)
\]
and the total PCG complexity is also
\[
O\!\big(N^d(\log N)^2\big)
\]
for \(d=1,2,3\) [2004.13961].

A distinct development concerns adaptive Legendre-Galerkin methods. In one dimension, the Babuška-Shen basis gives a coefficient representation of the \(H_0^1\)-norm and permits residual-based adaptivity for
\[
Lu := -D\cdot(\nu Du)+\sigma u = f \quad \text{in } I=(-1,1), \qquad u(-1)=u(1)=0.
\]
The Galerkin approximation on
\[
V_\Lambda=\operatorname{span}\{\eta_k:\,k\in\Lambda\}
\]
satisfies residual-error equivalence, and the ideal adaptive algorithm based on Dörfler marking yields the contraction
\[
|u-u_{n+1}|\le \rho\,|u-u_n|, \qquad
\rho(\theta)=\sqrt{1-\frac{\alpha_*}{\alpha^*}\theta^2}\in(0,1).
\]
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 [1206.5524].

In multiple dimensions, tensorized Babuška-Shen functions are no longer directly usable as an \(H_0^1(\Omega)\)-orthogonal basis. The multidimensional adaptive theory therefore introduces a quasi-orthonormalization that constructs a nearly-orthonormal Babuška-Shen basis \(\phi=\{\phi_k\}\) from the tensorized basis. This yields norm equivalences
\[
\|v\|_{H_0^1(\Omega)}^2 \simeq \sum_{k\in\mathcal K} d_k |\hat v_k|^2,
\qquad
\|f\|_{H^{-1}(\Omega)}^2 \simeq \sum_{k\in\mathcal K} d_k^{-1}|\hat f_k|^2,
\]
which enable feasible residual computation, Dörfler marking, enrichment, coarsening, and contraction:
\[
\|u-u_{n+1}\|_{H_0^1(\Omega)} \le \rho \|u-u_n\|_{H_0^1(\Omega)}.
\]
This places multidimensional adaptive Legendre-Galerkin discretization within the framework of nonlinear approximation and Gevrey-type sparsity [1408.0030].

## 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
\[
M\dot a=f(a),
\]
the Galerkin RONS modification enforces invariants \(I_k(a)\) by solving the constrained problem
\[
\min_{\dot a\in\mathbb R^N}\mathcal J(a,\dot a), \qquad \text{s.t.}\quad I_k(a(t))=\text{const.},
\]
which yields
\[
M\dot a = f(a)-\sum_{k=1}^m \lambda_k \nabla_a I_k(a),
\qquad
C(a)\lambda=b(a).
\]
The correction term is chosen so that
\[
\frac{d}{dt}I_k(a(t))=0
\]
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 \(\phi_i\) is Legendre or Legendre-like [2312.05322].

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,
\[
c_{N,T}\to 0,\qquad C_{N,T}\to\infty\qquad\text{as }N\to\infty.
\]
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 [1612.00332]. 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
\[
u_N(x)=\sum_{k=1}^{N-1}\alpha_k\,\phi_k(x),
\]
and the network is trained to predict \(\alpha\) by minimizing the residual of the discrete linear system
\[
(S+\nu M)\alpha=F.
\]
The population loss is exactly the squared algebraic residual, and the analysis proves that the empirical minimizer converges to the Legendre-Galerkin weak solution \(u_N\) under the stated assumptions [2211.08900]. 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 \(p\)-version SIPDG setting, the Legendre basis is modal, orthogonal in \(L^2\), 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 [2509.13669]. 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 [2211.10599], [2004.13961], [1612.00332].

Source: https://www.emergentmind.com/topics/legendre-galerkin-discretization