Elliptic PDE Solving Methods

Ground truth for the stationary finite-difference path. Elliptic problems have no time variable, so the boundary conditions are the problem: the solution is the direct solve of

\[ -\mathcal{T} u + \rho(x) u = f \]

rather than a backward march from terminal data. Here T is the transport operator built from Transport<N> and \(\rho\) is the reaction coefficient.

Relationship to the time-dependent contract

The parabolic and elliptic paths share the spatial discretization in finite_difference/discretization.rs. Both consume DimensionKind to select a stencil and Transport<N> to supply the per-dimension plus/minus coefficients and local source. They differ only in the outer solve strategy:

• PdeProblem<N> marches backward in time.
• EllipticProblem<N> assembles the operator once and solves it directly.

This keeps the stencils in one place while making the stationary strategy explicit. The stationary work lives inside finite_difference, not as a sibling module.

Elliptic problem contract

EllipticProblem<N> supplies:

• dimension_kind(dim): the stencil kind for coordinate dim.
• transport(state, derivs): the plus/minus transport coefficients. The returned source is ignored by the elliptic assembler; forcing belongs in rhs.
• rhs(state): the forcing f.
• reaction(state): a local zero-order coefficient added to the diagonal.
• boundary_conditions(): per-dimension lower and upper conditions.

The operator is assembled as

\[ -\sum_i \big[ a_i^{+} (u_{i,+} - u) + a_i^{-} (u_{i,-} - u) \big] + \rho(x)\, u = f(x), \]

where \(\rho\) is the reaction term. The negative sign follows the standard elliptic convention: for diffusion coefficient \(D\), \(T u = D u''\) and the problem reads \(-D u'' + \rho u = f\). The same DimensionKind stencil table as in solver_pde.md is used for \(a_i^{+}\) and \(a_i^{-}\).

Elliptic control problem contract

For infinite-horizon stochastic control, the solver consumes EllipticControlProblem<N>, which extends crate::models::control::ControlProblem<N> with the same transport /discretization data. Its stationary HJB equation is

\[ 0 = \sup_u \lbrace f(x,u) + \mathcal{L}^u V - r V \rbrace. \]

After optimizing over u, the equation is assembled in the same form as the pure elliptic problem:

\[ -\mathcal{T}^{u^*} V + r V = s^{u^*}, \]

where \(T V + s = \sup_u \lbrace f + \mathcal{L}^u V \rbrace\). The transport method returns plus, minus, and source = s; the discount \(r\) becomes the reaction diagonal.

StationarySolver::solve_control performs policy iteration:

1. Compute derivatives of the current value field.
2. Optimize the control at every node.
3. Assemble the operator and right-hand side from the optimized transport.
4. Solve the linear system.
5. Repeat until the value function stops changing.

The pure EllipticProblem path uses rhs for the forcing and ignores Transport::source; the control path uses Transport::source as the right-hand side.

Boundary conditions

BoundaryCondition is defined in discretization.rs:

• Dirichlet(value): \(u = \text{value}\).
• Neumann(value): \(du/dx = \text{value}\), with the derivative taken along the increasing coordinate direction \(x\), not an outward normal.
• Robin { value, alpha, beta }: \(\alpha u + \beta\, du/dx = \text{value}\).

BoundaryConditions<N> bundles the lower and upper conditions for every dimension. It defaults to Neumann(0) on every side.

Boundary condition implementation status

Solution strategies

StationarySolver assembles a sparse CSR operator and solves it:

• 1D problems use the tridiagonal LAPACK solve.
• N-dimensional problems use successive over-relaxation.

Dirichlet conditions replace the boundary row with the identity. Neumann and Robin conditions use a ghost point. The missing neighbour coefficient is folded into the existing interior neighbour and produces a right-hand-side (and diagonal, for Robin) correction. The derivative in Neumann and Robin is taken along the increasing coordinate direction \(x\).

Reference problems

The implementation is validated against closed forms in solver/tests/elliptic_benchmarks.rs (pure elliptic) and solver/tests/elliptic_control.rs (the control path). The exact solutions are stated once in Exact Solutions to PDEs; the problems are:

• 1D Laplace (Dirichlet), Poisson, and reaction-diffusion;
• 2D Poisson with separable forcing;
• 1D Laplace with Neumann and Robin lower boundaries.

These six problems cover constant, polynomial, separable multi-dimensional, reaction-dominated, Neumann, and Robin boundary behaviour, and verify that the assembler solves the bare elliptic operator \(-\mathcal{T} u + \rho u = f\) rather than an identity-shifted operator.

The control path is validated against the stationary (algebraic Riccati) LQ solution in Exact Solutions in Stochastic Optimal Control, with \(a = -1\), \(b = 1\), \(q = r = 1\), \(\rho = 0.1\), and \(c = 0\). Dirichlet boundary conditions match the exact value at the grid endpoints, so the test isolates the interior operator and the policy-iteration loop.