API reference

pdesolve(pde, grid, guess [, τs]; kwargs...)

Solve a system of nonlinear ODEs/PDEs — typically Hamilton–Jacobi–Bellman equations — by finite differences. Handles an arbitrary number of coupled unknown functions on a state grid with any positive number of state variables.

Positional arguments

• pde: the local equation, a function (state, u) -> out (or (state, u, t) -> out for a time-dependent equation), where
  • state is a NamedTuple with the current grid point (one entry per state variable),
  • u is a NamedTuple with each unknown function and its finite-difference derivatives at that point (e.g. v, vk_up, vk_down, vkk, …),
  • out is a NamedTuple with one time derivative per unknown (e.g. (; vt)).
  pde may also return a second NamedTuple of objects to save on the grid; these are returned in result.saved.
• grid: a NamedTuple (or OrderedDict) mapping each state variable name to an AbstractVector (its grid), e.g. (; k = range(...)).
• guess: a NamedTuple (or OrderedDict) mapping each unknown name to an array of initial values with the same shape as the grid. For a time-dependent problem, it is the terminal value at τs[end].
• τs (optional): an increasing time grid. The equation is then solved backward from τs[end], and each solution array gains a trailing time dimension: result.solution.v[.., i] is the solution at time τs[i].

Keyword arguments

• bc: boundary derivatives, a NamedTuple (or OrderedDict) mapping Symbol(unknown, state) to a (lower, upper) tuple (scalars, or arrays matching the boundary slice). Entries may be given for any subset of (unknown, state) pairs; omitted pairs default to reflecting boundaries (zero first derivative at the boundary). Both entries are ∂v/∂state oriented in the increasing-state direction.
• is_algebraic: a NamedTuple (or OrderedDict) of Bools marking equations that are algebraic (no time derivative) rather than PDEs. Defaults to all false.
• lower_bound, upper_bound: lower/upper bounds on the unknowns for HJB variational inequalities (optimal stopping), solved as a mixed complementarity problem. Pass either arrays in the same order as the unknowns or a NamedTuple/OrderedDict with the same names as guess. Default to -Inf/Inf (unbounded).
• alg = NonlinearSolve.NewtonRaphson(): NonlinearSolve algorithm used for each nonlinear solve. EconPDEs exports the NonlinearSolve module, so pass any compatible algorithm object, for example alg = NonlinearSolve.TrustRegion().
• abstol: convergence tolerance on the residual. Defaults to sqrt(eps()).
• maxiters: maximum number of pseudo-transient iterations. Defaults to 100.
• Δ: initial pseudo-transient time step for stationary problems. Defaults to 1.0; pass Δ = Inf to solve the stationary residual in one nonlinear solve, with no continuation. May also be a NamedTuple/OrderedDict giving one step per unknown (scalars or grid-sized arrays; omitted unknowns default to 1.0). The sign of an entry sets the direction of that unknown's false transient: positive marches its equation backward in time (the value-function convention vt = -(RHS - ρv)), negative forward — so a market-clearing equation written in its natural tâtonnement direction (rt = r_implied - r) relaxes toward equilibrium with a negative step, e.g. Δ = (; v = 1.0, r = -1.0). ±Inf entries disable damping for that unknown. The entries fix a relative profile; the continuation scales all of them by one common adaptive factor. Not accepted when a time grid τs is supplied; then the spacing of τs is the time step.
• verbose: print a one-line problem summary, convergence progress, and a final convergence summary. Defaults to true. With verbose = false a successful solve prints nothing — convergence failures are still reported with @warn — so pdesolve can run inside a loop (e.g. an estimation) without flooding the log.
• warn: report convergence failures with @warn even when verbose = false. Defaults to true. Pass warn = false when the caller checks result.converged itself and a failure is an expected outcome (e.g. probing whether a guess is already good enough for a one-shot Δ = Inf solve).
• Further stationary solver knobs (scale, minΔ, maxΔ, max_residual_growth) are forwarded to finiteschemesolve. Further inner-solve knobs (inner_maxiters, inner_abstol, inner_verbose) are also accepted for time-dependent problems.
• check_monotonicity: if true, warn when the assembled residual Jacobian has same-variable spatial off-diagonal entries with the wrong monotonicity sign. Defaults to false. A flipped relaxation direction — a residual-Jacobian diagonal opposite in sign to the unknown's pseudo-time step at every grid point, e.g. returning RHS - ρv instead of -(RHS - ρv) with a positive Δ — is detected per unknown and warned about regardless of this option (whenever the unknown's Δ entry is finite).
• monotonicity_tol: tolerance for the monotonicity check. Defaults to 1e-6.
• monotonicity_max_warnings: maximum number of monotonicity warnings to print. Defaults to 5.

Returns an EconPDEResult with fields solution (a NamedTuple with the solved unknowns), residual_norm, and saved (a NamedTuple with the saved objects, together with the solved unknowns, or an empty NamedTuple if the PDE saves nothing). result.converged reports whether the residual met the tolerance (across all times, for a time-dependent problem).

EconPDEResult

The value returned by pdesolve, with fields:

• solution: the solved unknown functions, a NamedTuple of arrays on the state grid. For a time-dependent problem each array has a trailing time dimension, so result.solution.v[.., i] is the solution at time τs[i].
• residual_norm: the norm of the residual at the solution (one per time step for a time-dependent problem).
• saved: objects saved by the PDE function, together with the solved unknowns, as a NamedTuple of arrays with the same layout as solution. Empty if the PDE saves nothing.

result.converged reports whether the residual norm is below the solver tolerance (the maximum across times, for a time-dependent problem).

finiteschemesolve(G!, y0; kwargs...)

Lower-level nonlinear solver behind pdesolve. Finds y such that G!(ydot, y) writes the residual into ydot and returns zero, using a nonlinear solver (Newton by default) with pseudo-transient continuation from the initial guess y0.

pdesolve assembles G! (the finite-difference residual of the PDE) and calls this function, so most users should call pdesolve instead. It returns the tuple (y, residual_norm).

Keyword arguments

• is_algebraic: Bool per entry of y0, marking algebraic (no time-derivative) equations.
• lower_bound, upper_bound: lower/upper bounds on y. When any bound is finite, the problem is solved with the minmax mixed-complementarity reformulation.
• abstol = sqrt(eps()): convergence tolerance on the residual norm.
• verbose = true: print one line per pseudo-transient time step (Iter, TimeStep, Residual) and a final convergence summary. A ✗ in place of the residual means the inner nonlinear solve failed at that Δ: no step is taken (so there is no residual to show) and it is retried with Δ/10 — the unusual causes, a NaN from the model or a singular Jacobian, are flagged in parentheses. With verbose = false a successful solve prints nothing; convergence failures are still reported with @warn.
• warn = true: report convergence failures with @warn even when verbose = false, so solves embedded in silenced loops (e.g. estimation) still surface failures. Pass warn = false when the caller inspects convergence itself and a failure is an expected outcome (e.g. probing whether a guess is good enough for a one-shot solve).
• alg = NonlinearSolve.NewtonRaphson(): NonlinearSolve algorithm used for each nonlinear solve. EconPDEs exports the NonlinearSolve module, so pass any compatible algorithm object, for example alg = NonlinearSolve.TrustRegion().
• maxiters = 100: maximum number of pseudo-transient (outer) iterations.
• Δ = 1.0: initial pseudo-transient time step. Δ = Inf solves the stationary residual in one nonlinear solve, with no continuation. May also be a vector with one entry per element of y0. The sign of an entry sets the direction of that equation's false transient: positive marches it backward in time (the value-function convention, which requires a positive residual-Jacobian diagonal), negative forward (e.g. a market-clearing equation written in its natural tâtonnement direction rt = r_implied - r, whose diagonal is negative); ±Inf entries disable damping for that equation. The entries fix a relative profile: the continuation scales all of them by one common adaptive factor.
• scale = 10.0: growth factor for the time step. After a successful step that reduces the residual, Δ is multiplied by scale * (old residual / new residual), and scale compounds across consecutive improving steps. After a failed inner solve, Δ is divided by 10 and its regrowth is capped at half the failed Δ, so a Δ that just failed is not immediately retried; the cap relaxes by scale with each subsequent successful step.
• minΔ = 1e-9, maxΔ = Inf: bounds on the time step. The solve stops when Δ < minΔ (typically a sign that the scheme is non-monotone or the initial guess is poor). For a vector Δ, the bounds — like the TimeStep column of the verbose output — refer to the smallest finite |Δ| entry.
• max_residual_growth = 10.0: guard on accepted steps. A step whose inner solve converges but whose stationary residual exceeds max_residual_growth times the previous residual (a blowup, or a NaN, which counts as infinite growth) is reverted and retried with a smaller Δ instead of accepted. The default is deliberately loose: a false transient may legitimately increase the residual along the way, and rejecting every worsening step can deadlock the continuation at minΔ.
• inner_maxiters = 10, inner_abstol = sqrt(eps()), inner_verbose = false: iteration limit, tolerance, and verbosity for each implicit time step (the inner solve).
• jac = nothing: optional in-place Jacobian of the stationary residual G!, with signature jac(J, y). EconPDEs adds the pseudo-time diagonal implied by Δ and is_algebraic before handing the Jacobian to NonlinearSolve.
• jac_prototype = nothing: prototype/sparsity pattern for the residual Jacobian. This is metadata for the NonlinearFunction, not an algorithm choice. When jac is omitted, a sparse prototype makes EconPDEs compute the Jacobian by colored finite differences; when jac is supplied, the prototype controls the allocated Jacobian type.
• colorvec = nothing: column coloring of jac_prototype for the colored finite differences (no two columns of the same color may share a nonzero row). Defaults to a greedy coloring of jac_prototype; pass one when it is known in closed form, as pdesolve does for its stencil pattern.