Edge, corner and vertex: Discretizing boundary conditions

Tianhao Zhao

Notation

In economic applications, we often deal with the following boundary problem of a nonlinear parabolic equation:

where $\mathcal{L}$$\mathcal{L}$ is the infinitesimal generator operator; $\mathcal{X}$$\mathcal{X}$ represents the interior of the hyper-rectangular domain; $\partial \mathcal{X}$$\partial\mathcal{X}$ denotes the boundary of the domain; $g_k(x)\mapsto \mathbb{R}$$g_k(x) \mapsto \mathbb{R}$ are evaluable functions; and Neumann boundary conditions are assumed, which align with the majority of our economic applications. In the previous blog post, we discussed the monotonicity of the numerical scheme, with a focus on the interior points of the domain. However, special attention must be paid when handling boundary points to ensure monotonicity, consistency and stability in the numerical solution.

Conventionally, we denote the discretized node space of finite points as ${\hat{\mathcal{X}}}_N\subset \mathcal{X}$$\hat{\mathcal{X}}_N \subset \mathcal{X}$, where $N$$N$ is the number of supporting nodes, and thus the boundary of such node space as $\partial {\hat{\mathcal{X}}}_N\subset \partial \mathcal{X}$$\partial\hat{\mathcal{X}}_N \subset \partial \mathcal{X}$. Usually, there are two strategies to handle the boundary $\partial \mathcal{X}$$\partial \mathcal{X}$ in the finite difference framework by discretizing the nodes in $\partial {\hat{\mathcal{X}}}_N$$\partial\hat{\mathcal{X}}_N$:

1. True boundary: Assuming all $x\in \partial {\hat{\mathcal{X}}}_N$$x\in\partial \hat{\mathcal{X}}_N$ are exactly on the boudary such that boundary conditions $v_k(x)=g_k(x)$$v_k(x) = g_k(x)$ hold for all these $x$$x$ nodes.

2. Ghost boundary: Assuming all $x\in \partial {\hat{\mathcal{X}}}_N$$x\in\partial\hat{\mathcal{X}}_N$ are slightly inward away from the boundary. Assuming there exist a "thin" suface $\partial {\hat{\mathcal{X}}}_N^g$$\partial\hat{\mathcal{X}}^g_N$ just outside $\partial {\hat{\mathcal{X}}}_N$$\partial\hat{\mathcal{X}}_N$ where the boundary conditions hold on points in $\partial {\hat{\mathcal{X}}}_N^g$$\partial\hat{\mathcal{X}}^g_N$ rather than on points in $\partial {\hat{\mathcal{X}}}_N$$\partial\hat{\mathcal{X}}_N$.

We discuss each of them while considering the quasi-implicit scheme on uniform/dense grid.

True boundary

Let $x[i]\in \partial {\hat{\mathcal{X}}}_N$$x[i] \in \partial\hat{\mathcal{X}}_N$ be a node that we are about to discretize, where $i$$i$ is a multi-index of $D$$D$ dimensions. We use $x_k[i]$$x_k[i]$ to denote the value of the $k$$k$-th dimension of $x[i]$$x[i]$. Before moving to the "corner" (intersection of two boundaries) or "vertex" (intersection of three or more boundaries), let's first discuss the "edge" case. Suppose $x[i]$$x[i]$ is on the boundary of only one dimension, say $m$$m$. For all $k\ne m$$k \neq m$, $x_k[i]$$x_k[i]$ are interior.

Edge

In the case of true boundary, if the boundary condition is Dirichlet, then we simply evaluate $v(x[i])$$v(x[i])$. However, for Neumann boundary conditions, we still need to discretize the boundary condition at $x[i]$$x[i]$. Recall the boundary condition along dimension $m$$m$:

The quasi-implicit method requires discretizing Eq ($\text{4}$$\ref{eq:04}$) on the layer of $t+{\mathrm{\Delta }}_t$$t+\Delta_t$. In this case, if $x[i]$$x[i]$ is the left boundary along dimension $m$$m$ i.e. $x_m[i]={\underset{―}{x}}^m$$x_m[i]=\underline{x}^m$, then we do forward difference and rearrange the terms to make all coefficients positive:

where the prime mark denotes the time layer of $t+{\mathrm{\Delta }}_t$$t +\Delta_t$. This scheme is unconditional monotonic. Similarly, in the case of right boundary, there is:

Equation ($\text{5},\text{6}$$\ref{eq:05},\ref{eq:06}$) are the equations for updating the guess of $v(x[i])$$v(x[i])$. One should NOT use the parabolic equation ($\text{1}$$\ref{eq:01}$) which holds on interior points.

Corner and vertex

In the case of the intersection of two boundaries, say $m,q=1,\dots ,D$$m, q = 1, \dots, D$ and $m\ne q$$m \neq q$, let $x[i]$$x[i]$ lie on the boundary along dimensions $m$$m$ and $q$$q$ simultaneously. A node like this is often referred to as a corner node. If three or more dimensions intersect simultaneously, we call such a node a vertex node. Note that the literature often does not distinguish between the two, which is fine.

The challenge in discretization here is to ensure that $v(x[i])$$v(x[i])$ simultaneously satisfies both boundary conditions. If we use Dirichlet conditions for both $m$$m$ and $q$$q$, we can always assign a consistent value to $v^{\prime }(x[i])$$v'(x[i])$. However, if two Neumann conditions meet, a problem arises: we have two equations but only one unknown, $v^{\prime }(x[i])$$v'(x[i])$. This indeterminacy persists even if there is no mixed partial derivative $v_{mq}(x[i])$$v_{mq}(x[i])$ involved.

One possible solution is to average the two Neumann conditions such that (assuming both boundaries are left boundaries):

The monotonicity still holds. This can trivially extends to three-dimension vertex nodes even nodes of higher dimensionality. However, one issue here is that the case of simultaneous-holding is certainly a solution to Eq. ($\text{7}$$\ref{eq:10}$), but it is obviously not the unique solution. Numerical non-monotonicity, instability and other violations may still exist. Therefore, it is critical to ensure that the implied value function $v(x[i])$$v(x[i])$ derived from the two Neumann boundary conditions is continuous at the corner node. In the context of a household problem, for example, consider the reflection boundary condition (Neumann) for liquid assets, which implies consuming all flow income, alongside the reflection boundary condition for illiquid housing wealth, which implies purchasing as much housing wealth as possible. These two conditions must determine an allocation at the “lowest consumption and lowest housing wealth” point: how should the flow income be allocated between consumption and housing wealth investment? Furthermore, is the manually assigned allocation consistent with the policy functions for consumption and housing wealth investment? Answering these questions requires a deep understanding of the behavior of the quantitative model before attempting computational implementation. In such cases, solving a discrete-time version of the problem on a coarse grid can provide valuable intuition. The discrete-time version and its continuous-time counterpart must exhibit consistency, while the discrete-time model typically does not require explicitly specified boundary conditions, instead relying on numerical optimization to determine solutions.

The idea of re-introducing some numerical optimization reminds us that, instead of manually specifying the allocation at such corner nodes, one can solve the optimization problem at $x[i]$$x[i]$ numerically. This approach preserves the convenience of analytical policy solutions for continuous-time models, as the number of corner/vertex nodes is very small relative to the total number of nodes. An example of such an optimization step looks like:

where $c$$c$ is a vector of control variables. Here, we are searching for a control that satisfies the maximum principle while respecting the multiple boundary conditions as much as possible.

It is important to note that, in economic applications, we typically do not impose "hard" boundary conditions that require strict equalities. Instead, we often work with "soft" boundary conditions, which are expressed as inequalities. The nature of the free boundary problem we are addressing validates the optimization framework outlined above without compromising adherence to the boundary conditions.

Example: optimization at corner nodes

Let's consider a simple human capital problem which has two states: $k$$k$ physical capital, and $h$$h$ human capital.

in which households choose consumption and the time allocation between working and education. We use soft reflection boundaries for $k$$k$ such that $\dot{k}(k,h)\ge 0$$\dot{k}(k,h) \geq 0$ for $k=\underset{―}{k}$$k = \underline{k}$, $\dot{k}(k,h)\le 0$$\dot{k}(k,h) \leq 0$ for $k=\overline{k}$$k = \bar{k}$, $\dot{h}(k,h)\ge 0$$\dot{h}(k,h) \geq 0$ for $h=\underset{―}{h}$$h = \underline{h}$, and $\dot{h}(k,h)\le 0$$\dot{h}(k,h) \leq 0$ for $h=\overline{h}$$h = \bar{h}$. Readers may notice that we use all "soft" boundary conditions rather than standard Neumann conditions because, by definition, this is a free-boundary problem.

Now, let's define the optimization problem. We define the generalized Hamiltonian $\mathcal{H}$$\mathcal{H}$, which is discounted to the current moment rather than to $t=0$$t = 0$ as in the standard statement of maximum principle.

Notes: The solved policy function of $c$$c$ and $s$$s$, which are denoted as functions of partial derivetives of $v$$v$, are:

$$\begin{array}{}\text{(17)}& & c=v_k^{-1}\text{(18)}& & s={\left(\frac{v_h}{v_k}\frac{\theta \alpha }{w}\right)}^{\frac{1}{1-\alpha }}\frac{1}{h}\end{array}$$

They will be discretized in Eq ($\text{24},\text{25}$$\ref{eq:24},\ref{eq:25}$). We list them here for reference.

For the "bottom-left" corner $x[i]$$x[i]$ in which $k=\underset{―}{k}$$k=\underline{k}$ and $h=\underset{―}{h}$$h=\underline{h}$, we discretize the Hamiltonian using forward difference along both dimensions (recall that we are handling the case of true boundary). Assume the values of neighbor nodes are known such that they become parameters of the Hamiltonian. The function value $v(x[i])$$v(x[i])$ becomes a control that we want to solve:

where $u:=v(x[i])$$u:=v(x[i])$ for convenience; $u^k:=v[i+{\mathrm{\Delta }}_k{e}_k]$$u^k := v[i+\Delta_k e_k]$ and $u^h:=v[i+{\mathrm{\Delta }}_h{e}_h]$$u^h:=v[i+\Delta_h e_h]$ are the two neighbor nodes used in the finite difference approximation. This discretized Hamiltonian, which serves as the objective function, is subject to:

where the two monotonicity constraints are imposed by the non-negativity constraint of consumption and our domain knowledge, respectively. Meanwhile, we require policy functions of $c$$c$ and $s$$s$ still hold in order to smooth the policy function surface, while make the optimal $u$$u$ bounded.

In this way, by knowing the two neighboring nodes $u^k$$u^k$ and $u^h$$u^h$, one can solve the above optimization problem efficiently using fmincon() or other nonlinear optimizer routines. This approach focuses solely on one point, $v(x[i])$$v(x[i])$, of the value function, allowing the problem to be modeled as a closure. Consequently, there is no need to interpolate the entire $v(\cdot )$$v(\cdot)$. Similarly, one can define such an optimization problem for all four corners in this human capital example. Furthermore, the problem can be trivially generalized to higher dimensions.

Remark: explicit? or implicit?

The optimization solution for corner/vertex nodes described above requires knowing the values of $v(x[i\pm {\mathrm{\Delta }}_k{e}_k])$$v(x[i \pm \Delta_k e_k])$ and $v(x[i\pm {\mathrm{\Delta }}_h{e}_h])$$v(x[i \pm \Delta_h e_h])$, i.e., the value function in the neighborhood of the corner nodes. This approach works perfectly with an explicit scheme as Eq ($\text{26}$$\ref{eq:23.5}$) since these values are guesses from the previous iteration. However, with the implicit method, these neighboring node values are on the time layer $t+{\mathrm{\Delta }}_t$$t + \Delta_t$, meaning they are unknown when solving for $v^{\prime }(x[i])$$v'(x[i])$. Equivalently, this results in nonlinear equation(s) that must be solved to update the corner nodes. Consequently, one needs to solve these nonlinear corner equations simultaneously with the interior equations. This process breaks the linearity of the quasi-implicit scheme that is commonly used.

However, one may have noticed that the above optimization problems are computationally very cheap to solve, making them suitable to be wrapped as a callable object. The corner/vertex equations can then be expressed as:

The monotonicity of this nonlinear scheme is always satisfied due to the monotonicity constraint embedded within the optimization routine. Meanwhile, the monotonicity of the scheme for interior nodes can be ensured by the regular quasi-implicit upwind scheme.

In general, either the implicit form as Eq ($\text{27}$$\ref{eq:23}$) or the explicit form as Eq ($\text{26}$$\ref{eq:23.5}$) works, regardless of which form is applied to interior nodes. Mixing implicit and explicit methods for different parts of the domain $\mathcal{X}$$\mathcal{X}$ is feasible and is often referred to as a locally defined difference scheme.

In practice, the explicit form is usually preferred due to its computational simplicity. However, if the implicit form is chosen, the entire numerical scheme (incorporating both interior and boundary nodes) remains implicit. To solve the entire system of updates, one may consider employing iterative algorithms such as the Newton method.

Hint: Always check CFL condition carefully when working with any explicit scheme.

Ghost boundary

The ghost boundary method is more frequently used due to its smoothness in the neighborhood of $\partial \mathcal{X}$$\partial \mathcal{X}$, especially when we are not particularly interested in the boundary behavior of the PDE but merely use the boundary to pin down the solution.

Remark: In Benjamin Moll's illustrative examples, the partial derivatives from backward difference at a left boundary node are manually specified (e.g., using the reflection of zero savings). This approach effectively embeds a backward difference operation with a ghost outside node.

The idea of discretizing ghost boundaries is similar to the true boundary case. However, in this approach, we require the PDE to hold on $x[i]\in \partial {\hat{\mathcal{X}}}_N$$x[i] \in \partial\hat{\mathcal{X}}_N$ while the boundary conditions are enforced on ghost nodes, e.g., on $x[i-{\mathrm{\Delta }}_k{e}_k]\in \partial {\hat{\mathcal{X}}}_N^g$$x[i - \Delta_k e_k] \in \partial\hat{\mathcal{X}}^g_N$.

Edge

Let's illustrate this using an edge node $x[i]$$x[i]$ which is the left boundary of dimension $h$$h$. Let's consider Neumann boundary condition on ghost node $v_m(x[i-{\mathrm{\Delta }}_m{e}_m])=g_m(x[i-{\mathrm{\Delta }}_m{e}_m])$$v_m(x[i-\Delta_me_m])=g_m(x[i-\Delta_me_m])$. Because it is a left boundary, we use forward difference to approximate $v_m(x[i-{\mathrm{\Delta }}_m{e}_m])$$v_m(x[i-\Delta_me_m])$, rearrange terms and get the following formula for $v(x[i-{\mathrm{\Delta }}_m{e}_m])$$v(x[i-\Delta_m e_m])$:

Suppose the PDE for interior points bas the following simple form:

Discretizing it and apply quasi-implcit upwind scheme:

It is easy to realize Eq ($\text{32}$$\ref{eq:28}$) is unconditionally monotonic.

Corner and vertex

When using ghost boundaries, a notable advantage is that we usually do not need to solve for the very corner ghost node unless we are approximating specific partial derivatives. To illustrate, consider $x[i]$$x[i]$, a "bottom-left" corner node. The following table provides the spatial locations of $x[i]$$x[i]$ and its neighboring ghost nodes:

The ghost corner, $x[i-{\mathrm{\Delta }}_x{e}_x-{\mathrm{\Delta }}_y{e}_y]$$x[i - \Delta_x e_x - \Delta_y e_y]$, is only required when approximating the off-diagonal element of the Hessian matrix, such as $v_{xy}(x[i])$$v_{xy}(x[i])$, particularly using backward difference. This type of approximation typically arises, at least in the context of economic research, in portfolio selection models with systemic risks. The same concept can be trivially generalized to vertex nodes.

However, if someone really wants to solve $v(x[i-{\mathrm{\Delta }}_x{e}_x-{\mathrm{\Delta }}_y{e}_y])$$v(x[i - \Delta_x e_x - \Delta_y e_y])$, they may simply use the optimization idea introduced in the last section. Notably, this optimization does not contribute to updating the guess in this example.