Method Overview

This page gives a compact user-facing overview of the rectangular-domain method implemented by gbmsc-pde. It is intended to explain the objects used by the package, not to replace the full mathematical reference. For citation guidance, see citation.md.

The current implementation targets rectangular domains, structured tensor-product source grids, and linear PDE rows evaluated at source nodes.

1. Grid and Local Subgrids

Let

\[ \Omega = [a,b] \times [c,d] \]

and let the structured grid be

\[ \mathcal{X}_h = \left\{ (x_i,y_j) : 0 \le i < N_x,\; 0 \le j < N_y \right\}. \]

For a uniform grid,

\[ x_i = a + i h_x, \qquad y_j = c + j h_y, \]

where

\[ h_x = \frac{b-a}{N_x-1}, \qquad h_y = \frac{d-c}{N_y-1}. \]

The method builds local tensor-product subgrids. Each local subgrid has shape (n_x, n_y) and contains the source nodes used by one local tensor-product interpolant. The step=(step_x, step_y) parameter controls how these subgrids are sampled across the source grid. When the step is smaller than the subgrid size, neighboring subgrids overlap; a source node can then belong to several nearby subgrids. This overlap is what allows the local interpolants to be blended smoothly by the Shepard weights. The set of sampled subgrids is denoted by

\[ \mathcal{G}_h. \]

2. Local Tensor-Product Lagrange Interpolation

For each sampled subgrid \(s \in \mathcal{G}_h\), let

\[ x^{(s)}_0,\ldots,x^{(s)}_{n_x-1}, \qquad y^{(s)}_0,\ldots,y^{(s)}_{n_y-1} \]

be its one-dimensional node coordinates. The local Lagrange basis functions are

\[ \ell^{(s)}_\alpha(x), \qquad \eta^{(s)}_\beta(y). \]

For source values \(u_{\alpha\beta}^{(s)}\) on the subgrid, the local tensor-product interpolant is

\[ I_s u(x,y) = \sum_{\alpha=0}^{n_x-1} \sum_{\beta=0}^{n_y-1} u_{\alpha\beta}^{(s)} \ell^{(s)}_\alpha(x) \eta^{(s)}_\beta(y). \]

The first derivatives are

\[ \partial_x I_s u(x,y) = \sum_{\alpha,\beta} u_{\alpha\beta}^{(s)} \frac{d\ell^{(s)}_\alpha}{dx}(x) \eta^{(s)}_\beta(y), \]

and

\[ \partial_y I_s u(x,y) = \sum_{\alpha,\beta} u_{\alpha\beta}^{(s)} \ell^{(s)}_\alpha(x) \frac{d\eta^{(s)}_\beta}{dy}(y). \]

The second derivatives are computed similarly using second derivatives of the one-dimensional Lagrange bases.

3. Shepard Blending

Each local interpolant is blended by a Shepard weight. Let \(z=(x,y)\) and let \(z_k^{(s)}\) be the source nodes in subgrid \(s\). With Shepard exponent \(\mu > 2\), the raw subgrid weight is represented by

\[ q_s(z) = \prod_{z_k^{(s)} \in s} \|z-z_k^{(s)}\|^{-\mu}. \]

The original normalized Shepard weight uses all sampled subgrids in the denominator:

\[ D(z) = \sum_{r\in\mathcal{G}_h} q_r(z), \qquad p_s(z) = \frac{q_s(z)}{D(z)}. \]

The blended Grid-Based Multinode Shepard Collocation approximation is

\[ U_h(z) = \sum_{s\in\mathcal{G}_h} p_s(z) I_su(z). \]

For interpolation at arbitrary query points, the package uses all sampled subgrids:

\[ \mathcal{S}(z) = \mathcal{G}_h. \]

For PDE rows evaluated at source nodes, the package also provides three nodal support modes through GBMSCApproximation. support_mode="active" uses the sampled subgrids that contain the source node and applies the original epsilon-regularized nodal formula. support_mode="all" keeps active contributors in the numerator but normalizes with all sampled subgrids. support_mode="nodal_limit" evaluates the source-node limiting formula by cancelling the common singular factor shared by active subgrids; this is the recommended mode for new rectangular PDE experiments. See gbs_approximation.md for usage details.

4. Nodal Differential Operators

The rectangular LinearPDE solver assembles sparse matrices acting on the flattened source vector

\[ \mathbf{u} = \left(u(z_0),u(z_1),\ldots,u(z_{N-1})\right)^T. \]

The derivative matrices are the sparse nodal operators used by the PDE and boundary-condition assembly:

• \(M_x\) maps nodal values to an approximation of \(\partial_x u\) at the source nodes.
• \(M_y\) maps nodal values to an approximation of \(\partial_y u\) at the source nodes.
• \(M_{xx}\) maps nodal values to an approximation of \(\partial_{xx}u\) at the source nodes.
• \(M_{yy}\) maps nodal values to an approximation of \(\partial_{yy}u\) at the source nodes.

Equivalently, these matrices satisfy

\[ M_x \mathbf{u} \approx \partial_x u \quad \text{at source nodes}, \]

\[ M_y \mathbf{u} \approx \partial_y u \quad \text{at source nodes}, \]

and

\[ M_{xx}\mathbf{u} \approx \partial_{xx}u, \qquad M_{yy}\mathbf{u} \approx \partial_{yy}u. \]

Because the LinearPDE rows are evaluated at source nodes, the tensor-product Lagrange basis has a nodal structure. The assembled derivative rows are sparse: each row only touches source nodes in local subgrids associated with that row node.

5. Boundary Conditions

Boundary conditions are registered on flattened rectangular grid-node ids.

For Dirichlet data,

\[ u(z_i) = g_i, \]

the corresponding matrix row is replaced by an identity row.

For Neumann data,

\[ \frac{\partial u}{\partial n}(z_i) = \nabla u(z_i)\cdot n_i = g_i, \]

the row is assembled as

\[ n_{x,i} M_x[i,:] + n_{y,i} M_y[i,:]. \]

For Robin data,

\[ \alpha_i u(z_i) + \beta_i \frac{\partial u}{\partial n}(z_i) = g_i, \]

the row is

\[ \alpha_i e_i^T + \beta_i \left( n_{x,i} M_x[i,:] + n_{y,i} M_y[i,:] \right), \]

where \(e_i\) is the unit coordinate row associated with node \(i\).