Solve nonlinear PDEs arising from economic models
Julia
143
469 commits
updated Aug 28, 2026
EconPDEs.jl solves nonlinear ODEs and PDEs that arise in economic models, especially Hamilton-Jacobi-Bellman equations.
You specify three things: (i) the state space, (ii) an initial guess for the solution, and (iii) a function that returns the PDE residual at each grid point. The package handles the rest — it builds the finite-difference derivatives, applies upwind stencils, assembles the sparse Jacobian, and solves the system by pseudo-transient Newton iteration.
Use it for stationary or time-dependent HJBs on one or more state variables, with one or more value functions.
EconPDEs.jl is robust, fast, and solves virtually all textbook continuous-time models; see the Examples for a list.
The package is registered in the Julia General registry:
] add EconPDEs
Current versions require Julia 1.10 or later.
Here is a small linear equation on one state, with a mean-reverting drift:
\rho v(x) = x + \mu(x) v_x(x),
\qquad
\mu(x) = -\kappa (x - 0.5).
pdesolve finds the stationary solution by embedding the equation in a time-dependent
problem, solved backward in time from an initial guess until the time derivative is zero.
\rho v(x, t) = x + \mu(x) v_x(x, t) + v_t(x, t),
Here is how to solve this equation using this package:
using EconPDEs
const rho = 0.05
const kappa = 0.2
# 1. Define a `NamedTuple` representing the state space as a grid:
grid = (; x = range(0.0, 1.0, length = 50))
# 2. Define a `NamedTuple` representing an initial guess for the solution:
guess = (; v = zeros(length(grid.x)))
# 3. Define a function encoding the PDE at each grid point:
function pde(state::NamedTuple, u::NamedTuple)
# Its first argument holds the coordinates of the current grid point.
(; x) = state
# Its second argument holds the function and its (finite-difference) derivatives at that point.
(; v, vx_up, vx_down, vxx) = u
# The function must return the time derivative implied by rho*v = x + mu*v_x + v_t. Note we upwind the first derivative.
mu = -kappa * (x - 0.5)
vt = -(x + mu * (mu >= 0 ? vx_up : vx_down) - rho * v)
return (; vt)::NamedTuple
end
# 4. Pass these three objects to `pdesolve`
result = pdesolve(pde, grid, guess)
# Solving for 1 unknown (v) on a 50-point grid
# Iter TimeStep Residual
# ---- -------- --------
# 1 1.00e+00 7.51e-02
# 2 1.09e+01 4.45e-02
# 3 1.85e+03 4.67e-04
# 4 1.76e+08 1.22e-10
# Converged after 4 time steps (63ms)
# EconPDEResult
# solution: v (50)
# residual_norm: 1.22e-10
# converged: true (tolerance 1.49e-08)
Note that, in encoding the pde, we used the forward difference when the drift is positive, the backward one otherwise (upwinding).
The same interface extends to multiple states...
grid = (; x = range(0.0, 1.0, length = 50),
y = collect(range(0.5, 1.5, length = 10)))
guess = (; v = zeros(length(grid.x), length(grid.y)))
function pde(state, u)
(; x, y) = state
(; v, vx_up, vx_down, vy_up, vy_down, vxx, vyy, vxy_up, vxy_down) = u
...
return (; vt)
end
... as well as multiple value functions (i.e., coupled PDES).
grid = (; x = range(0.0, 1.0, length = 50))
guess = (; v = zeros(length(grid.x)), w = zeros(length(grid.x)))
function pde(state, u)
(; x) = state
(; v, vx_up, vx_down, vxx, w, wx_up, wx_down, wxx) = u
...
return (; vt, wt)
end
Read the documentation for further details.
The documentation also includes a gallery of runnable examples (see examples/ for the original scripts):
Please cite EconPDEs.jl if it is useful in your work. You can use the "Cite this repository" button on the GitHub page, or the following BibTeX entry:
@software{EconPDEs.jl,
author = {Gomez, Matthieu},
title = {{EconPDEs.jl}: {Solving PDEs in economics}},
url = {https://github.com/matthieugomez/EconPDEs.jl},
}
Julia
100.0%
Solve nonlinear PDEs arising from economic models
Julia
143
469 commits
updated Aug 28, 2026
EconPDEs.jl solves nonlinear ODEs and PDEs that arise in economic models, especially Hamilton-Jacobi-Bellman equations.
You specify three things: (i) the state space, (ii) an initial guess for the solution, and (iii) a function that returns the PDE residual at each grid point. The package handles the rest — it builds the finite-difference derivatives, applies upwind stencils, assembles the sparse Jacobian, and solves the system by pseudo-transient Newton iteration.
Use it for stationary or time-dependent HJBs on one or more state variables, with one or more value functions.
EconPDEs.jl is robust, fast, and solves virtually all textbook continuous-time models; see the Examples for a list.
The package is registered in the Julia General registry:
] add EconPDEs
Current versions require Julia 1.10 or later.
Here is a small linear equation on one state, with a mean-reverting drift:
\rho v(x) = x + \mu(x) v_x(x),
\qquad
\mu(x) = -\kappa (x - 0.5).
pdesolve finds the stationary solution by embedding the equation in a time-dependent
problem, solved backward in time from an initial guess until the time derivative is zero.
\rho v(x, t) = x + \mu(x) v_x(x, t) + v_t(x, t),
Here is how to solve this equation using this package:
using EconPDEs
const rho = 0.05
const kappa = 0.2
# 1. Define a `NamedTuple` representing the state space as a grid:
grid = (; x = range(0.0, 1.0, length = 50))
# 2. Define a `NamedTuple` representing an initial guess for the solution:
guess = (; v = zeros(length(grid.x)))
# 3. Define a function encoding the PDE at each grid point:
function pde(state::NamedTuple, u::NamedTuple)
# Its first argument holds the coordinates of the current grid point.
(; x) = state
# Its second argument holds the function and its (finite-difference) derivatives at that point.
(; v, vx_up, vx_down, vxx) = u
# The function must return the time derivative implied by rho*v = x + mu*v_x + v_t. Note we upwind the first derivative.
mu = -kappa * (x - 0.5)
vt = -(x + mu * (mu >= 0 ? vx_up : vx_down) - rho * v)
return (; vt)::NamedTuple
end
# 4. Pass these three objects to `pdesolve`
result = pdesolve(pde, grid, guess)
# Solving for 1 unknown (v) on a 50-point grid
# Iter TimeStep Residual
# ---- -------- --------
# 1 1.00e+00 7.51e-02
# 2 1.09e+01 4.45e-02
# 3 1.85e+03 4.67e-04
# 4 1.76e+08 1.22e-10
# Converged after 4 time steps (63ms)
# EconPDEResult
# solution: v (50)
# residual_norm: 1.22e-10
# converged: true (tolerance 1.49e-08)
Note that, in encoding the pde, we used the forward difference when the drift is positive, the backward one otherwise (upwinding).
The same interface extends to multiple states...
grid = (; x = range(0.0, 1.0, length = 50),
y = collect(range(0.5, 1.5, length = 10)))
guess = (; v = zeros(length(grid.x), length(grid.y)))
function pde(state, u)
(; x, y) = state
(; v, vx_up, vx_down, vy_up, vy_down, vxx, vyy, vxy_up, vxy_down) = u
...
return (; vt)
end
... as well as multiple value functions (i.e., coupled PDES).
grid = (; x = range(0.0, 1.0, length = 50))
guess = (; v = zeros(length(grid.x)), w = zeros(length(grid.x)))
function pde(state, u)
(; x) = state
(; v, vx_up, vx_down, vxx, w, wx_up, wx_down, wxx) = u
...
return (; vt, wt)
end
Read the documentation for further details.
The documentation also includes a gallery of runnable examples (see examples/ for the original scripts):
Please cite EconPDEs.jl if it is useful in your work. You can use the "Cite this repository" button on the GitHub page, or the following BibTeX entry:
@software{EconPDEs.jl,
author = {Gomez, Matthieu},
title = {{EconPDEs.jl}: {Solving PDEs in economics}},
url = {https://github.com/matthieugomez/EconPDEs.jl},
}
Julia
100.0%