The Multivariate Theory of Connections ^†

Abstract

This paper extends the univariate Theory of Connections, introduced in (Mortari, 2017), to the multivariate case on rectangular domains with detailed attention to the bivariate case. In particular, it generalizes the bivariate Coons surface, introduced by (Coons, 1984), by providing analytical expressions, called constrained expressions, representing all possible surfaces with assigned boundary constraints in terms of functions and arbitrary-order derivatives. In two dimensions, these expressions, which contain a freely chosen function, $g(x,y)$, satisfy all constraints no matter what the $g(x,y)$ is. The boundary constraints considered in this article are Dirichlet, Neumann, and any combinations of them. Although the focus of this article is on two-dimensional spaces, the final section introduces the Multivariate Theory of Connections, validated by mathematical proof. This represents the multivariate extension of the Theory of Connections subject to arbitrary-order derivative constraints in rectangular domains. The main task of this paper is to provide an analytical procedure to obtain constrained expressions in any space that can be used to transform constrained problems into unconstrained problems. This theory is proposed mainly to better solve PDE and stochastic differential equations.

1. Introduction

The Theory of Connections (ToC), as introduced in [1], consists of a general analytical framework to obtain constrained expressions, $f\left(x\right)$, in one-dimension. A constrained expression is a function expressed in terms of another function, $g\left(x\right)$, that is freely chosen and, no matter what the $g\left(x\right)$ is, the resulting expression always satisfies a set of n constraints. ToC generalizes the one-dimensional interpolation problem subject to n constraints using the general form,

$${\displaystyle f\left(x\right)=g\left(x\right)+\sum _{k=1}^n{\eta }_k{p}_k\left(x\right),}$$

(1)

where $p_k\left(x\right)$ are n user-selected linearly independent functions, ${\eta }_k$ are derived by imposing the n constraints, and $g\left(x\right)$ is a freely chosen function subject to be defined and nonsingular where the constraints are specified. Besides this requirement, $g\left(x\right)$ can be any function, including, discontinuous functions, delta functions, and even functions that are undefined in some domains. Once the ${\eta }_k$ coefficients have been derived, then Equation (1) satisfies all the n constraints, no matter what the $g\left(x\right)$ function is.

Constrained expressions in the form given in Equation (1) are provided for a wide class of constraints, including constraints on points and derivatives, linear combinations of constraints, as well as infinite and integral constraints [2]. In addition, weighted constraints [3] and point constraints on continuous and discontinuous periodic functions with assigned period can also be obtained [1]. How to extend ToC to inequality and nonlinear constraints is currently a work in progress.

The Theory of Connections framework can be considered the generalization of interpolation; rather than providing a class of functions (e.g., monomials) satisfying a set of n constraints, it derives all possible functions satisfying the n constraints by spanning all possible $g\left(x\right)$ functions. This has been proved in Ref. [1]. A simple example of a constrained expression is,

$$f\left(x\right)=g\left(x\right)+{\displaystyle \frac{x(2x_2-x)}{2(x_2-x_1)}}\left[{\dot{y}}_1-\dot{g}\left(x_1\right)\right]+{\displaystyle \frac{x(x-2x_1)}{2(x_2-x_1)}}\left[{\dot{y}}_2-\dot{g}\left(x_2\right)\right].$$

(2)

This equation always satisfies ${\left.{\displaystyle \frac{\mathrm{d}f}{\mathrm{d}x}}\right|}_{x_1}={\dot{y}}_1$ and ${\left.{\displaystyle \frac{\mathrm{d}f}{\mathrm{d}x}}\right|}_{x_2}={\dot{y}}_2$, as long as $\dot{g}\left(x_1\right)$ and $\dot{g}\left(x_2\right)$ are defined and nonsingular. In other words, the constraints are embedded into the constrained expression.

Constrained expressions can be used to transform constrained optimization problems into unconstrained optimization problems. Using this approach, fast least-squares solutions of linear [4] and nonlinear [5] ODE have been obtained at machine error accuracy and with low (actually, very low) condition number. Direct comparisons of ToC versus MATLAB’s ode45 [6] and Chebfun [7] have been performed on a small test of ODE with excellent results [4,5]. In particular, the ToC approach to solve ODE consists of a unified framework to solve IVP, BVP, and multi-value problems. The extension of differential equations subject to component constraints [8] has opened the possibility for ToC to solve in real-time a class of direct optimal control problems [9], where the constraints connect state and costate.

This study first extends the Theory of Connections to two-dimensions by providing, for rectangular domains, all surfaces that are subject to: (1) Dirichlet constraints; (2) Neumann constraints; and (3) any combination of Dirichlet and Neumann constraints. This theory is then generalized to the Multivariate Theory of Connections which provide in n-dimensional space all possible manifolds that satisfy boundary constraints on the value and boundary constraints on any-order derivative.

This article is structured as follows. First, it shows that the one-dimensional ToC can be used in two dimensions when the constraints (functions or derivatives) are provided along one axis only. This is a particular case, where the original univariate theory [1] can be applied with basically no modifications. Then, a two dimensional ToC version is developed for Dirichlet type boundary constraints. This theory is then extended to include Neumann and mixed type boundary constraints. Finally, the theory is extended to n-dimensions and to incorporate arbitrary-order derivative boundary constraints followed by a mathematical proof validating it.

2. Manifold Constraints in One Axis, Only

Consider the function, $f\left(\mathit{x}\right)$, where $f:{\mathbb{R}}^n\to {\mathbb{R}}^1$, subject to one constraint manifold along the ith variable, $x_i$, that is, ${f\left(\mathit{x}\right)|}_{x_i=v}=c\left({\mathit{x}}_i^v\right)$. For instance, in 3-D space, this can be the surface constraint, ${f(x,y,z)|}_{y=\pi }=c(x,\pi ,z)$. All manifolds satisfying this constraint can be expressed using the additive form provided in Ref. [1],

$$f\left(\mathit{x}\right)=g\left(\mathit{x}\right)+\left[c\left({\mathit{x}}_i^v\right)-g\left({\mathit{x}}_i^v\right)\right]$$

where $g\left(\mathit{x}\right)$ is a freely chosen function that must be defined and nonsingular at the constraint coordinates. When m manifold constraints are defined along the $x_i$-axis, then the 1-D methodology [1] can be applied as it is. For instance, the constrained expression subject to m constraints along the $x_i$ variable evaluated at $x_i=w_k$, where $k\in [1,m]$, that is, ${f\left(\mathit{x}\right)|}_{x_i=w_k}=c\left({\mathit{x}}_i^{w_k}\right)$, is,

$${\displaystyle f\left(\mathit{x}\right)=g\left(\mathit{x}\right)+\sum _{k=1}^m\left\{{\displaystyle \left[c\left({\mathit{x}}_i^{w_k}\right)-g\left({\mathit{x}}_i^{w_k}\right)\right]\prod _{j\ne k}{\displaystyle \frac{x_i-w_j}{w_k-w_j}}}\right\}.}$$

(3)

Note that this equation coincides with the Waring interpolation form (better known as Lagrangian interpolation form) [10] if the free function vanishes, $g\left(\mathit{x}\right)=0$.

2.1. Example #1: Surface Subject to Four Function Constraints

The first example is designed to show how to use Equation (3) with mixed, continuous, discontinuous, and multiple constraints. Consider the following four constraints,

$$c(x,-2)=sin\left(2x\right),c(x,0)=3cosx\left[(x+1)\mathrm{mod}\left(2\right)\right],c(x,1)=9e^{-x^2},\mathrm{and}c(x,3)=1-x.$$

This example highlights that the constraints and free-function may be discontinuous by using the modular arithmetic function. The result is a surface that is continuous in x at some coordinates (at y = −2, 1, and 3) and discontinuous at $y=0$. The surfaces shown in Figure 1 and Figure 2 were obtained using two distinct expressions for the free function, $g(x,y)$.

2.2. Example #2: Surface Subject to Two Functions and One Derivative Constraint

This second example is provided to show how to use the general approach given in Equation (1) and described in [1], when derivative constraints are involved. Consider the following three constraints,

$$c(x,-2)=sin\left(2x\right),c_y(x,0)=0,\mathrm{and}c(x,1)=9e^{-x^2}.$$

Using the functions $p_1\left(y\right)=1$, $p_2\left(y\right)=y$, and $p_3\left(y\right)=y^2$, the constrained expression form satisfying these three constraints assumes the form,

$$f(x,y)=g(x,y)+{\eta }_1\left(x\right)+{\eta }_2\left(x\right)y+{\eta }_3\left(x\right)y^2.$$

(4)

The three constraints imply the constraints,

$$\begin{array}{ccc}\hfill sin\left(2x\right)& =& g(x,-2)+{\eta }_1-2{\eta }_2+4{\eta }_3\hfill \\ \hfill 0& =& g_y(x,0)+{\eta }_2\hfill \\ \hfill 9e^{-x^2}& =& g(x,1)+{\eta }_1+{\eta }_2+{\eta }_3,\hfill \end{array}$$

from which the values of the ${\eta }_k$ coefficients,

$$\begin{array}{ccc}\hfill {\eta }_1& =& 2g_y(x,0)+12e^{-x^2}-{\displaystyle \frac{sin\left(2x\right)}{3}}+{\displaystyle \frac{1}{3}}g(x,-2)-{\displaystyle \frac{4}{3}}g(x,1)\hfill \\ \hfill {\eta }_2& =& -g_y(x,0)\hfill \\ \hfill {\eta }_3& =& {\displaystyle \frac{sin\left(2x\right)}{3}}-{\displaystyle \frac{1}{3}}g(x,-2)-g_y(x,0)-3e^{-x^2}+{\displaystyle \frac{1}{3}}g(x,1),\hfill \end{array}$$

can be derived. After substituting these coefficients into Equation (4), the constrained expression that always satisfies the three initial constraints is obtained. Using this expression and two different free functions, $g(x,y)$, we obtained the surfaces shown in Figure 3 and Figure 4, respectively. The constraint $c_y(x,0)=0$, difficult to see in both figures, can be verified analytically.

3. Connecting Functions in Two Directions

In this section, the Theory of Connections is extended to the two-dimensional case. Note that dealing with constraints in two (or more) directions (functions or derivatives) requires particular attention. In fact, two orthogonal constraint functions cannot be completely distinct as they intersect at one point where they need to match in value. In addition, if the formalism derived for the 1-D case is applied to 2-D case, some complications arise. These complications are highlighted in the following simple clarifying example.

Consider the two boundary constraint functions, $f(x,0)=q\left(x\right)$ and $f(0,y)=h\left(y\right)$. Searching the constrained expression as originally done for the one-dimensional case implies the expression,

$$f(x,y)=g(x,y)+{\eta }_1{p}_1(x,y)+{\eta }_2{p}_2(x,y).$$

The constraints imply the two constraints,

$$\left\{\begin{array}{c}q\left(x\right)=g(x,0)+{\eta }_1{p}_1(x,0)+{\eta }_2{p}_2(x,0)\hfill \\ h\left(y\right)=g(0,y)+{\eta }_1{p}_1(0,y)+{\eta }_2{p}_2(0,y).\hfill \end{array}\right.$$

To obtain the values of ${\eta }_1$ and ${\eta }_2$, the determinant of the matrix to invert is $p_1(x,0)p_2(0,y)-p_1(0,y)p_2(x,0)$. This determinant is y by selecting $p_1(x,y)=1$ and $p_2(x,y)=y$, or it is x by selecting $p_1(x,y)=x$ and $p_2(x,y)=1$. Therefore, to avoid singularities, this approach requires paying particular attention to the domain definition and/or on the user-selected functions, $p_k(x,y)$. To avoid dealing with these issues, a new (equivalent) formalism to derive constrained expressions is devised for the higher dimensional case.

The Theory of Connections extension to the higher dimensional case (with constraints on all axes) can be obtained by re-writing the constrained expression into an equivalent form, highlighting a general and interesting property. Let us show this by an example. Equation (2) can be re-written as,

$$f\left(x\right)=\underset{A\left(x\right)}{\underbrace{{\displaystyle \frac{x(2x_2-x)}{2(x_2-x_1)}}{\dot{y}}_1+{\displaystyle \frac{x(x-2x_1)}{2(x_2-x_1)}}{\dot{y}}_2}}+\underset{B\left(x\right)}{\underbrace{g\left(x\right)-{\displaystyle \frac{x(2x_2-x)}{2(x_2-x_1)}}{\dot{g}}_1-{\displaystyle \frac{x(x-2x_1)}{2(x_2-x_1)}}{\dot{g}}_2}}.$$

(5)

These two components, $A\left(x\right)$ and $B\left(x\right)$, of a constrained expression have a specific general meaning. The term, $A\left(x\right)$, represents an (any) interpolating function satisfying the constraints while the $B\left(x\right)$ term represents all interpolating functions that are vanishing at the constraints. Therefore, the generation of all functions satisfying multiple orthogonal constraints in n-dimensional space can always be expressed by the general form, $f\left(\mathit{x}\right)=A\left(\mathit{x}\right)+B\left(\mathit{x}\right)$, where $A\left(\mathit{x}\right)$ is any function satisfying the constraints and $B\left(\mathit{x}\right)$ must represent all functions vanishing at the constraints. Equation $f\left(\mathit{x}\right)=A\left(\mathit{x}\right)+B\left(\mathit{x}\right)$ is actually an alternative general form to write a constrained expression, that is, an alternative way to generalize interpolation: rather than derive a class of functions (e.g., monomials) satisfying a set of constraints, it represents all possible functions satisfying the set of constraints.

To prove that this additive formalism can describe all possible functions satisfying the constraints is immediate. Let $f\left(\mathit{x}\right)$ be all functions satisfying the constraints and $y\left(\mathit{x}\right)=A\left(\mathit{x}\right)+B\left(\mathit{x}\right)$ be the sum of a specific function satisfying the constraints, $A\left(\mathit{x}\right)$, and a function, $B\left(\mathit{x}\right)$, representing all functions that are null at the constraints. Then, $y\left(\mathit{x}\right)$ will be equal to $f\left(\mathit{x}\right)$ iff $B\left(\mathit{x}\right)=f\left(\mathit{x}\right)-A\left(\mathit{x}\right)$, representing all functions that are null at the constraints.

As shown in Equation (5), once the $A\left(\mathit{x}\right)$ function is obtained, then the $B\left(\mathit{x}\right)$ function can be immediately derived. In fact, $B\left(\mathit{x}\right)$ can be obtained by subtracting the $A\left(\mathit{x}\right)$ function, where all the constraints are specified in terms of the $g\left(\mathit{x}\right)$ free function, from the free function $g\left(\mathit{x}\right)$. For this reason, let us write the general expression of a constrained expression as,

$$f\left(\mathit{x}\right)=A\left(\mathit{x}\right)+g\left(\mathit{x}\right)-A\left(g\right(\mathit{x}\left)\right),$$

(6)

where $A\left(g\right(\mathit{x}\left)\right)$ indicates the function satisfying the constraints where the constraints are specified in term of $g\left(\mathit{x}\right)$.

The previous discussion serves to prove that the problem of extending Theory of Connections to higher dimensional spaces consists of the problem of finding the function, $A\left(\mathit{x}\right)$, only. In two dimensions, the function $A\left(\mathit{x}\right)$ is provided in literature by the Coons surface [11], $f(x,y)$. This surface satisfies the Dirichlet boundary constraints,

$$f(0,y)=c(0,y),f(1,y)=c(1,y),f(x,0)=c(x,0),\mathrm{and}f(x,1)=c(x,1),$$

(7)

where the surface is contained in the $x,y\in [0,1]\times [0,1]$ domain. This surface is used in computer graphics and in computational mechanics applications to smoothly join other surfaces together, particularly in finite element method and boundary element method, to mesh problem domains into elements. The expression of the Coons surface is,

$$\begin{array}{cc}\hfill f(x,y)=& (1-x)c(0,y)+xc(1,y)+(1-y)c(x,0)+yc(x,1)-xyc(1,1)\hfill \\ \hfill & -(1-x)(1-y)c(0,0)-(1-x)yc(0,1)-x(1-y)c(1,0),\hfill \end{array}$$

where the four subtracting terms are there for continuity. Note the constraint functions at boundary corners must have the same value, $c(0,0)$, $c(0,1)$, $c(1,0)$, and $c(1,1)$. This equation can be written in matrix form as,

$$f(x,y)=\left\{\begin{array}{ccc}1,& 1-x,& x\end{array}\right\}\left[\begin{array}{ccc}\hfill 0& c(x,0)& c(x,1)\hfill \\ \hfill c(0,y)& -c(0,0)& -c(0,1)\hfill \\ \hfill c(1,y)& -c(1,0)& -c(1,1)\hfill \end{array}\right]\left\{\begin{array}{c}1\\ 1-y\\ y\end{array}\right\},$$

or, equivalently,

$$f(x,y)={\mathit{v}}^T\left(x\right)\mathcal{M}\left(c(x,y)\right)\mathit{v}\left(y\right),$$

(8)

where

$$\mathcal{M}\left(c(x,y)\right)=\left[\begin{array}{ccc}\hfill 0& c(x,0)& c(x,1)\hfill \\ \hfill c(0,y)& -c(0,0)& -c(0,1)\hfill \\ \hfill c(1,y)& -c(1,0)& -c(1,1)\hfill \end{array}\right]\mathrm{and}\mathit{v}\left(z\right)=\left\{\begin{array}{c}1\\ 1-z\\ z\end{array}\right\}.$$

Since the $f(x,y)$ boundaries match the boundaries of the $c(x,y)$ constraint function, then the identity, $f(x,y)={\mathit{v}}^T\left(x\right)\mathcal{M}\left(f(x,y)\right)\mathit{v}\left(y\right)$, holds for any $f(x,y)$ function. Therefore, the $B\left(\mathit{x}\right)$ function can be set as,

$$B\left(\mathit{x}\right):=g(x,y)-{\mathit{v}}^T\left(x\right)\mathcal{M}\left(g(x,y)\right)\mathit{v}\left(y\right),$$

(9)

representing all functions that are always zero at the boundary constraints, as $g(x,y)$ is a free function.

4. Theory of Connections Surface Subject to Dirichlet Constraints

Equations (8) and (9) can be merged to provide all surfaces with the boundary constraints defined in Equation (7) in the following compact form,

$$f(x,y)=\underset{A(x,y)}{\underbrace{{\mathit{v}}^T\left(x\right)\mathcal{M}\left(c(x,y)\right)\mathit{v}\left(y\right)}}+\underset{B\left(\mathit{x},y\right)}{\underbrace{g(x,y)-{\mathit{v}}^T\left(x\right)\mathcal{M}\left(g(x,y)\right)\mathit{v}\left(y\right)}}.$$

(10)

where, again, $A(x,y)$ indicates an expression satisfying the boundary function constraints defined by $c(x,y)$ and $B(x,y)$ an expression that is zero at the boundaries. In matrix form, Equation (10) becomes,

$$f(x,y)={\left\{\begin{array}{c}1\\ 1-x\\ x\end{array}\right\}}^T\left[\begin{array}{ccc}g(x,y)& c(x,0)-g(x,0)& c(x,1)-g(x,1)\\ c(0,y)-g(0,y)& g(0,0)-c(0,0)& g(0,1)-c(0,1)\\ c(1,y)-g(1,y)& g(1,0)-c(1,0)& g(1,1)-c(1,1)\end{array}\right]\left\{\begin{array}{c}1\\ 1-y\\ y\end{array}\right\},$$

where $g(x,y)$ is a freely chosen function. In particular, if $g(x,y)=0$, then the ToC surface becomes the Coons surface.

Figure 5 (left) shows the Coons surface subject to the constraints,

$$\begin{array}{cc}\hfill c(x,0)& =sin(3x-\pi /4)cos(\pi /3)\hfill \\ \hfill c(x,1)& =sin(3x-\pi /4)cos(4+\pi /3)\hfill \\ \hfill c(0,y)& =sin(-\pi /4)cos(4y+\pi /3)\hfill \\ \hfill c(1,y)& =sin(3-\pi /4)cos(4y+\pi /3),\hfill \end{array}$$

and Figure 5 (right) shows a ToC surface that is obtained using the free function,

$$g(x,y)={\displaystyle \frac{1}{3}}cos\left(4\pi x\right)sin\left(6\pi y\right)-x^{2}cos\left(2\pi y\right).$$

(11)

For generic boundaries defined in the rectangle $x,y\in [x_i,x_f]\times [y_i,y_f]$, the ToC surface becomes,

$$\begin{array}{cc}\hfill f(x,y)=& g(x,y)+{\displaystyle \frac{x-x_f}{x_i-x_f}}\left[c(x_i,y)-g(x_i,y)\right]+{\displaystyle \frac{x-x_i}{x_f-x_i}}\left[c(x_f,y)-g(x_f,y)\right]\hfill \\ \hfill & +{\displaystyle \frac{y-y_f}{y_i-y_f}}\left[c(x,y_i)-g(x,y_i)\right]+{\displaystyle \frac{y-y_i}{y_f-y_i}}\left[c(x,y_f)-g(x,y_f)\right]\hfill \\ \hfill & -{\displaystyle \frac{(x-x_f)(y-y_f)}{(x_i-x_f)(y_i-y_f)}}\left[c(x_i,y_i)-g(x_i,y_i)\right]\hfill \\ \hfill & -{\displaystyle \frac{(x-x_f)(y-y_i)}{(x_i-x_f)(y_f-y_i)}}\left[c(x_i,y_f)-g(x_i,y_f)\right]\hfill \\ \hfill & -{\displaystyle \frac{(x-x_i)(y-y_f)}{(x_f-x_i)(y_i-y_f)}}\left[c(x_f,y_i)-g(x_f,y_i)\right]\hfill \\ \hfill & -{\displaystyle \frac{(x-x_i)(y-y_i)}{(x_f-x_i)(y_f-y_i)}}\left[c(x_f,y_f)-g(x_f,y_f)\right].\hfill \end{array}$$

(12)

Equation (12) can also be set in matrix form,

$$f(x,y)={\mathit{v}}_x^T(x,x_i,x_f)\mathcal{M}(x,y){\mathit{v}}_y(y,y_i,y_f)$$

where

$$\mathcal{M}(x,y)=\left[\begin{array}{ccc}g(x,y)& c(x,y_i)-g(x,y_i)& c(x,y_f)-g(x,y_f)\\ c(x_i,y)-g(x_i,y)& g(x_i,y_i)-c(x_i,y_i)& g(x_i,y_f)-c(x_i,y_f)\\ c(x_f,y)-g(x_f,y)& g(x_f,y_i)-c(x_f,y_i)& g(x_f,y_f)-c(x_f,y_f)\end{array}\right]$$

and

$${\mathit{v}}_x(x,x_i,x_f)=\left\{\begin{array}{c}1\\ {\displaystyle \frac{x-x_f}{x_i-x_f}}\\ {\displaystyle \frac{x-x_i}{x_f-x_i}}\end{array}\right\}\mathrm{and}{\mathit{v}}_y(y,y_i,y_f)=\left\{\begin{array}{c}1\\ {\displaystyle \frac{y-y_f}{y_i-y_f}}\\ {\displaystyle \frac{y-y_i}{y_f-y_i}}\end{array}\right\}.$$

Note that all the ToC surfaces provided are linear in $g(x,y)$, and, therefore, they can be used to solve, by linear/nonlinear least-squares, two-dimensional optimization problems subject to boundary function constraints, such as linear/nonlinear partial differential equations.

5. Multi-Function Constraints at Generic Coordinates

Equation (12) can be generalized to many function constraints (grid of functions). Assume a set of $n_x$ function constraints $c(x_k,y)$ and a set of $n_y$ function constraints $c(x,y_k)$ intersecting at the $n_x{n}_y$ points $p_{ij}=c(x_i,y_j)$, then all surfaces satisfying the $n_x{n}_y$ function constraints can be expressed by,

$$\begin{array}{cc}\hfill f(x,y)=g(x,y)& {\displaystyle +\sum _{k=1}^{n_x}\left[c(x_k,y)-g(x_k,y)\right]\prod _{i\ne k}{\displaystyle \frac{x-x_i}{x_k-x_i}}}\hfill \\ \hfill & {\displaystyle +\sum _{k=1}^{n_y}\left[c(x,y_k)-g(x,y_k)\right]\prod _{i\ne k}{\displaystyle \frac{y-y_i}{y_k-y_i}}}\hfill \\ \hfill & {\displaystyle -\sum _{i=1}^{n_x}\left\{{\displaystyle \sum _{j=1}^{n_y}{\displaystyle \frac{(x-x_j)(y-y_i)}{(x_i-x_j)(y_j-y_i)}}\left[c(x_i,y_j)-g(x_i,y_j)\right]}\right\}.}\hfill \end{array}$$

(13)

Again, Equation (13) can be written in compact form,

$$f(x,y)={\mathit{v}}^T\left(x\right)\mathcal{M}\left(c(x,y)\right)\mathit{v}\left(y\right)+g(x,y)-{\mathit{v}}^T\left(x\right)\mathcal{M}\left(g(x,y)\right)\mathit{v}\left(y\right)$$

where,

$$\mathit{v}\left(x\right)=\left\{\begin{array}{c}1\\ {\displaystyle \prod _{i\ne 1}{\displaystyle \frac{x-x_i}{x_1-x_i}}}\\ ⋮\\ {\displaystyle \prod _{i\ne n_x}{\displaystyle \frac{x-x_i}{x_{n_x}-x_i}}}\end{array}\right\}\mathrm{and}\mathit{v}\left(y\right)=\left\{\begin{array}{c}1\\ {\displaystyle \prod _{i\ne 1}{\displaystyle \frac{y-y_i}{y_1-y_i}}}\\ ⋮\\ {\displaystyle \prod _{i\ne n_y}{\displaystyle \frac{y-y_i}{y_{n_y}-y_i}}}\end{array}\right\}$$

and

$$\mathcal{M}\left(c(x,y)\right)=\left[\begin{array}{cccc}\hfill 0& c(x,y_1)& \dots \hfill & c(x,y_{n_y})\\ \hfill c(x_1,y)& -c(x_1,y_1)& \dots \hfill & -c(x_1,y_{N_y})\\ \hfill ⋮& ⋮& \ddots \hfill & ⋮\\ \hfill c(x_{n_x},y)& -c(x_{n_x},y_1)& \dots \hfill & -c(x_{n_x},y_{n_y})\end{array}\right]$$

For example, two function constraints in x and three function constraints in y can be obtained using the matrix,

$$\mathcal{M}\left(c(x,y)\right)=\left[\begin{array}{cccc}\hfill 0& c(x,y_1)& c(x,y_2)\hfill & c(x,y_3)\\ \hfill c(x_1,y)& -c(x_1,y_1)& -c(x_1,y_2)\hfill & -c(x_1,y_3)\\ \hfill c(x_2,y)& -c(x_2,y_1)& -c(x_2,y_2)\hfill & -c(x_2,y_3)\end{array}\right]$$

and the vectors,

$$\mathit{v}\left(x\right)=\left\{\begin{array}{c}1\\ {\displaystyle \frac{x-x_2}{x_1-x_2}}\\ {\displaystyle \frac{x-x_1}{x_2-x_1}}\end{array}\right\}\mathrm{and}\mathit{v}\left(y\right)=\left\{\begin{array}{c}1\\ {\displaystyle \frac{(y-y_2)(y-y_3)}{(y_1-y_2)(y_1-y_3)}}\\ {\displaystyle \frac{(y-y_1)(y-y_3)}{(y_2-y_1)(y_2-y_3)}}\\ {\displaystyle \frac{(y-y_2)(y-y_1)}{(y_3-y_2)(y_3-y_1)}}\end{array}\right\}.$$

Two examples of ToC surfaces are given in Figure 6 in the $x,y\in [-2,1]\times [1,3]$ domain.

6. Constraints on Function and Derivatives

The “Boolean sum formulation” was provided by Farin [12] (also called “Hermite–Coons formulation”) of the Coons surface that includes boundary derivatives,

$$f(x,y)={\mathit{v}}^T\left(y\right){\mathit{F}}^x\left(x\right)+{\mathit{v}}^T\left(x\right){\mathit{F}}^y\left(y\right)-{\mathit{v}}^T\left(x\right)M^{xy}\mathit{v}\left(y\right)$$

(14)

where

$$\begin{array}{cc}\hfill \mathit{v}\left(z\right)& :={\{2z^3-3z^2+1,z^3-2z^2+z,-2z^3+3z^2,z^3-z^2\}}^T\hfill \\ \hfill {\mathit{F}}^x\left(x\right)& :={\{c(x,0),c_y(x,0),c(x,1),c_y(x,1)\}}^T\hfill \\ \hfill {\mathit{F}}^y\left(y\right)& :={\{c(0,y),c_x(0,y),c(1,y),c_x(1,y)\}}^T\hfill \end{array}$$

and

$$M^{xy}(x,y):=\left[\begin{array}{cccc}\hfill c(0,0)& c_y(0,0)& c(0,1)\hfill & c_y(0,1)\\ \hfill c_x(0,0)& c_{xy}(0,0)& c_x(0,1)\hfill & c_{xy}(0,1)\\ \hfill c(1,0)& c_y(1,0)& c(1,1)\hfill & c_y(1,1)\\ \hfill c_x(1,0)& c_{xy}(1,0)& c_x(1,1)\hfill & c_{xy}(1,1)\end{array}\right].$$

The formulation provided in Equation (14) can be put in the matrix compact form,

$$f(x,y)={\mathit{v}}^T\left(x\right)\mathcal{M}\left(c(x,y)\right)\mathit{v}\left(y\right),$$

(15)

where

$$\mathit{v}\left(z\right):={\{1,2z^3-3z^2+1,z^3-2z^2+z,-2z^3+3z^2,z^3-z^2\}}^T$$

(16)

and the $5\times 5$ matrix, $\mathcal{M}\left(c\right(x,y\left)\right)$, has the expression,

$$\mathcal{M}\left(c(x,y)\right):=\left[\begin{array}{ccccc}\hfill 0& c(x,0)& c_y(x,0)\hfill & c(x,1)& c_y(x,1)\\ \hfill c(0,y)& -c(0,0)& -c_y(0,0)\hfill & -c(0,1)& -c_y(0,1)\\ \hfill c_x(0,y)& -c_x(0,0)& -c_{xy}(0,0)\hfill & -c_x(0,1)& -c_{xy}(0,1)\\ \hfill c(1,y)& -c(1,0)& -c_y(1,0)\hfill & -c(1,1)& -c_y(1,1)\\ \hfill c_x(1,y)& -c_x(1,0)& -c_{xy}(1,0)\hfill & -c_x(1,1)& -c_{xy}(1,1)\end{array}\right].$$

(17)

To verify the boundary derivative constraints, the following partial derivatives of Equation (15) are used,

$$\begin{array}{cc}\hfill f_x(x,y)& =[{\mathit{v}}_x^T\left(x\right)\mathcal{M}\left(c(x,y)\right)+{\mathit{v}}^T\left(x\right){\mathcal{M}}_x\left(c(x,y)\right)]\mathit{v}\left(y\right)\hfill \\ \hfill f_y(x,y)& ={\mathit{v}}^T\left(x\right)[{\mathcal{M}}_y^T\left(c(x,y)\right)\mathit{v}\left(y\right)+\mathcal{M}\left(c(x,y)\right){\mathit{v}}_y\left(y\right)],\hfill \end{array}$$

where

$${\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}z}}=\left\{\begin{array}{c}0\\ 6z(z-1)\\ 3z^2-4z+1\\ 6z(1-z)\\ z(3z-2)\end{array}\right\},{\mathcal{M}}_y=\left[\begin{array}{cc}0& {\mathbf{0}}_{1\times 4}\\ c_y(0,y)& {\mathbf{0}}_{1\times 4}\\ c_{xy}(0,y)& {\mathbf{0}}_{1\times 4}\\ c_y(1,y)& {\mathbf{0}}_{1\times 4}\\ c_{xy}(1,y)& {\mathbf{0}}_{1\times 4}\end{array}\right],\mathrm{and}{\mathcal{M}}_x^T=\left[\begin{array}{cc}0& {\mathbf{0}}_{1\times 4}\\ c_x(x,0)& {\mathbf{0}}_{1\times 4}\\ c_{xy}(x,0)& {\mathbf{0}}_{1\times 4}\\ c_x(x,1)& {\mathbf{0}}_{1\times 4}\\ c_{xy}(x,1)& {\mathbf{0}}_{1\times 4}\end{array}\right].$$

The ToC in 2D with function and derivative boundary constraints is simply,

$$f(x,y)=\underset{A(x,y)}{\underbrace{{\mathit{v}}^T\left(x\right)\mathcal{M}\left(c(x,y)\right)\mathit{v}\left(y\right)}}+\underset{B(x,y)}{\underbrace{g(x,y)-{\mathit{v}}^T\left(x\right)\mathcal{M}\left(g(x,y)\right)\mathit{v}\left(y\right)}}$$

(18)

where the $\mathcal{M}$ matrix and the $\mathit{v}$ vectors are provided by Equations (17) and (16), respectively.

Dirichlet/Neumann mixed constraints can be derived, as shown in the examples provided in Section 6.1 through Section 6.4. The matrix compact form is simply obtained from the matrix defined in Equation (17) by removing the rows and the columns associated with the boundary constraints not provided, while the vectors $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ are derived by specifying the constraints. Note that in general the vectors $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ are not unique. The reason the vectors $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ are not unique comes from the fact that the $A\left(\mathit{x}\right)$ term in Equation (6) is not unique.

In the next subsections, four Dirichlet/Neumann mixed constraint examples providing the simplest expressions for $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ are derived. The Appendix A contains the expressions for the $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ vectors associated with all the combinations of Dirichlet and Neumann constraints.

6.1. Constraints: $c(0,y)$ and $c(x,0)$

In this case, the Coons-type surface satisfying the boundary constraints can be expressed as,

$$f(x,y)=\left\{\begin{array}{cc}1& p\left(x\right)\end{array}\right\}\left[\begin{array}{cc}\hfill 0& c(x,0)\\ \hfill c(0,y)& -c(0,0)\end{array}\right]\left\{\begin{array}{c}1\\ q\left(y\right)\end{array}\right\}$$

where $p\left(x\right)$ and $q\left(y\right)$ are unknown functions. Expanding, we obtain $f(x,y)=c(x,0)q\left(y\right)+p\left(x\right)\left[c\right(0,y)-c(0,0\left)q\right(y\left)\right]$. The two constraints are satisfied if,

$$\begin{array}{cc}\hfill c(0,y)& =c(0,0)q\left(y\right)+p\left(0\right)\left[c\right(0,y)-c(0,0\left)q\right(y\left)\right]\hfill \\ \hfill c(x,0)& =c(x,0)q\left(0\right)+p\left(x\right)\left[c\right(0,0)-c(0,0\left)q\right(0\left)\right].\hfill \end{array}$$

Therefore, the $p\left(x\right)$ and $q\left(y\right)$ functions must satisfy $p\left(0\right)=1$ and $q\left(0\right)=1$. The simplest expressions satisfying these equations can be obtained by selecting $p\left(x\right)=1$ and $q\left(y\right)=1$. In this case, the associated ToC surface is given by,

$$f(x,y)=\left\{\begin{array}{cc}1& 1\end{array}\right\}\left[\begin{array}{cc}g(x,y)& c(x,0)-g(x,0)\\ c(0,y)-g(0,y)& g(0,0)-c(0,0)\end{array}\right]\left\{\begin{array}{c}1\\ 1\end{array}\right\}$$

Note that any functions satisfying $p\left(0\right)=1$ and $q\left(0\right)=1$ can be adopted to obtain the ToC surface satisfying the constraints $f(0,y)=c(0,y)$ and $f(x,0)=c(x,0)$. This is because there are infinite Coons-type surfaces satisfying the constraints. Consequently, the vectors $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ are not unique.

6.2. Constraints: $c(0,y)$ and $c_y(x,0)$

For these boundary constraints, the Coons-type surface is expressed by,

$$f(x,y)=\left\{\begin{array}{cc}1& p\left(x\right)\end{array}\right\}\left[\begin{array}{cc}\hfill 0& \hfill c_y(x,0)\\ \hfill c(0,y)& -c_y(0,0)\hfill \end{array}\right]\left\{\begin{array}{c}1\\ q\left(y\right)\end{array}\right\}$$

= c_y(x, 0)q(y) + p(x)[c (0, y) − c_y (0, 0) q(y)].

The constraints are satisfied if,

$$\begin{array}{cc}\hfill c(0,y)& =c_y(0,0)q\left(y\right)+p\left(0\right)[c(0,y)-c_y(0,0)q\left(y\right)],\hfill \\ \hfill c_y(x,0)& =c_y(x,0)q_y\left(0\right)+p\left(x\right)[c_y(0,0)-c_y(0,0)q_y\left(0\right)].\hfill \end{array}$$

Therefore, the $p\left(x\right)$ and $q\left(y\right)$ functions must satisfy $p\left(0\right)=1$ and $q_y\left(0\right)=1$. One solution is $p\left(x\right)=1$ and $q\left(y\right)=y$. Therefore, the associated ToC surface is given by,

$$f(x,y)=\left\{\begin{array}{cc}1& 1\end{array}\right\}\left[\begin{array}{cc}g(x,y)& c_y(x,0)-g_y(x,0)\\ c(0,y)-g(0,y)& g_y(0,0)-c_y(0,0)\end{array}\right]\left\{\begin{array}{c}1\\ y\end{array}\right\}.$$

6.3. Neumann Constraints: $c_x(0,y)$, $c_x(1,y)$, $c_y(x,0)$, and $c_y(x,1)$

In this case, the Coons-type surface satisfying the boundary constraints can be expressed as,

$$f(x,y)=\left\{\begin{array}{ccc}1,& p_1\left(x\right),& p_2\left(x\right)\end{array}\right\}\left[\begin{array}{ccc}\hfill 0& c_y(x,0)& c_y(x,1)\hfill \\ \hfill c_x(0,y)& -c_{xy}(0,0)& -c_{xy}(0,1)\hfill \\ \hfill c_x(1,y)& -c_{xy}(1,0)& -c_{xy}(1,1)\hfill \end{array}\right]\left\{\begin{array}{c}1\\ q_1\left(y\right)\\ q_2\left(y\right)\end{array}\right\}.$$

The constraints are satisfied if,

$$\begin{array}{cc}\hfill c_x(0,y)& =q_1\left(y\right)c_{xy}(0,0)+q_2\left(y\right)c_{xy}(0,1)+\hfill \\ \hfill & +p_{1x}\left(0\right)[c_x(0,y)-q_1\left(y\right)c_{xy}(0,0)-q_2\left(y\right)c_{xy}(0,1)]+\hfill \\ \hfill & +p_{2x}\left(0\right)[c_x(1,y)-q_1\left(y\right)c_{xy}(1,0)-q_2\left(y\right)c_{xy}(1,1)]\hfill \\ \hfill c_x(1,y)& =q_1\left(y\right)c_{xy}(1,0)+q_2\left(y\right)c_{xy}(1,1)+\hfill \\ \hfill & +p_{1x}\left(1\right)[c_x(0,y)-q_1\left(y\right)c_{xy}(0,0)-q_2\left(y\right)c_{xy}(0,1)]+\hfill \\ \hfill & +p_{2x}\left(1\right)[c_x(1,y)-q_1\left(y\right)c_{xy}(1,0)-q_2\left(y\right)c_{xy}(1,1)]\hfill \\ \hfill c_y(x,0)& =q_{1y}\left(0\right)c_y(x,0)+q_{2y}\left(0\right)c_y(x,1)+\hfill \\ \hfill & +p_1\left(x\right)[c_{xy}(0,0)-q_{1y}\left(0\right)c_{xy}(0,0)-q_{2y}\left(0\right)c_{xy}(0,1)]+\hfill \\ \hfill & +p_2\left(x\right)[c_{xy}(1,0)-q_{1y}\left(0\right)c_{xy}(1,0)-q_{2y}\left(0\right)c_{xy}(1,1)]\hfill \\ \hfill c_y(x,1)& =q_{1y}\left(1\right)c_y(x,0)+q_{2y}\left(1\right)c_y(x,1)+\hfill \\ \hfill & +p_1\left(x\right)[c_{xy}(0,1)-q_{1y}\left(1\right)c_{xy}(0,0)-q_{2y}\left(1\right)c_{xy}(0,1)]+\hfill \\ \hfill & +p_2\left(x\right)[c_{xy}(1,1)-q_{1y}\left(1\right)c_{xy}(1,0)-q_{2y}\left(1\right)c_{xy}(1,1)].\hfill \end{array}$$

These equations imply $p_{1x}\left(0\right)=q_{1x}\left(0\right)=1$, $p_{1x}\left(1\right)=q_{1x}\left(1\right)=0$, $p_{2x}\left(0\right)=q_{2x}\left(0\right)=0$, and $p_{2x}\left(1\right)=q_{2x}\left(1\right)=1$. Therefore, the simplest solution is $p_1\left(t\right)=q_1\left(t\right)=t-t^2/2$ and $p_2\left(t\right)=q_2\left(t\right)=t^2/2$. Then, the associated ToC surface satisfying the Neumann constraints is given by,

$$f(x,y)={\mathit{v}}^T\left(x\right)\left[\begin{array}{ccc}g(x,y)& c_y(x,0)-g_y(x,0)& c_y(x,1)-g_y(x,1)\\ c_x(0,y)-g_x(0,y)& g_{xy}(0,0)-c_{xy}(0,0)& g_{xy}(0,1)-c_{xy}(0,1)\\ c_x(1,y)-g_x(1,y)& g_{xy}(1,0)-c_{xy}(1,0)& g_{xy}(1,1)-c_{xy}(1,1)\end{array}\right]\mathit{v}\left(y\right)$$

where

$${\mathit{v}}^T\left(x\right)=\left\{\begin{array}{ccc}1,& x-{\displaystyle \frac{x^2}{2}},& {\displaystyle \frac{x^2}{2}}\end{array}\right\}\mathrm{and}\mathit{v}\left(y\right)=\left\{\begin{array}{ccc}1,& y-{\displaystyle \frac{y^2}{2}},& {\displaystyle \frac{y^2}{2}}\end{array}\right\}.$$

6.4. Constraints: $c(0,y)$, $c_y(x,0)$, and $c_y(x,1)$

In this case, the Coons-type surface satisfying the boundary constraints is in the form,

$$f(x,y)={\left\{\begin{array}{c}1\\ p\left(x\right)\end{array}\right\}}^T\left[\begin{array}{ccc}\hfill 0& c_y(x,0)& c_y(x,1)\hfill \\ \hfill c(0,y)& -c_y(0,0)& -c_y(0,1)\hfill \end{array}\right]\left\{\begin{array}{c}1\\ q_1\left(y\right)\\ q_2\left(y\right)\end{array}\right\}.$$

The constraints are satisfied if $p\left(0\right)=1$, $p_{1y}\left(0\right)=1$, $p_{1y}\left(1\right)=0$, $p_{2y}\left(0\right)=0$, and $p_{2y}\left(1\right)=1$. Therefore, the associated ToC surface is,

$$f(x,y)={\left\{\begin{array}{c}1\\ 1\end{array}\right\}}^T\left[\begin{array}{ccc}g(x,y)& c_y(x,0)-g_y(x,0)& c_y(x,1)-g_y(x,1)\\ c(0,y)-g(0,y)& g_y(0,0)-c_y(0,0)& g_y(0,1)-c_y(0,1)\end{array}\right]\left\{\begin{array}{c}1\\ y-{\displaystyle \frac{y^2}{2}}\\ {\displaystyle \frac{y^2}{2}}\end{array}\right\}.$$

6.5. Generic Mixed Constraints

Consider the case of mixed constraints,

$$\begin{array}{c}f(x,y_1)=c(x,y_1)\\ f_x(x,y_2)=c_x(x,y_2)\\ f(x,y_3)=c(x,y_3)\end{array}\mathrm{and}\begin{array}{c}{f}_y(x_1,y)=c_y(x_1,y)\\ f_y(x_2,y)=c_y(x_2,y)\\ f(x_3,y)=c(x_3,y)\end{array}.$$

(19)

In this case, the surface satisfying the boundary constraints is built using the matrix,

$$\mathcal{M}\left(c(x,y)\right)=\left[\begin{array}{cccc}\hfill 0& c(x,y_1)& c_x(x,y_2)\hfill & c(x,y_3)\\ \hfill c_y(x_1,y)& -c_y(x_1,y_1)& -c_{xy}(x_1,y_2)\hfill & -c_y(x_1,y_3)\\ \hfill c_y(x_2,y)& -c_y(x_2,y_1)& -c_{xy}(x_2,y_2)\hfill & -c_y(x_2,y_3)\\ \hfill c(x_3,y)& -c(x_3,y_1)& -c_x(x_3,y_2)\hfill & -c(x_3,y_3)\end{array}\right]$$

and all surfaces subject to the constraints defined in Equation (19) can be obtained by,

$$f(x,y)=\mathit{v}{\left(x\right)}^T\mathcal{M}\left(c(x,y)\right)\mathit{v}\left(y\right)+g(x,y)-\mathit{v}{\left(x\right)}^T\mathcal{M}\left(g(x,y)\right)\mathit{v}\left(y\right),$$

where

$$\mathit{v}\left(x\right)=\left\{\begin{array}{c}1\\ p_1(x,x_1,x_2,x_3)\\ p_2(x,x_1,x_2,x_3)\\ p_3(x,x_1,x_2,x_3)\end{array}\right\}\mathrm{and}\mathit{v}\left(y\right)=\left\{\begin{array}{c}1\\ q_1(y,y_1,y_2,y_3)\\ q_2(y,y_1,y_2,y_3)\\ q_3(y,y_1,y_2,y_3)\end{array}\right\}$$

are vectors made of the (not unique) function vectors $\mathit{v}\left(x\right)$ and $\mathit{v}\left(y\right)$ whose expressions can be found by satisfying the constraints (as done in the previous four subsections) along with a methodology similar to that given in Section 5.

7. Extension to $\mathit{n}$-Dimensional Spaces and Arbitrary-Order Derivative Constraints

This section provides the Multivariate Theory of Connections, as the generalization to n-dimensional rectangular domains with arbitrary-order boundary derivatives of what is presented above for two-dimensional space. Using tensor notation, this generalization is represented in the following compact form,

$$F\left(\mathit{x}\right)=\underset{A\left(\mathit{x}\right)}{\underbrace{\mathcal{M}{\left(c\left(\mathit{x}\right)\right)}_{i_1{i}_2\dots i_n}{\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {\mathit{v}}_{i_n}}}+\underset{B\left(\mathit{x}\right)}{\underbrace{g\left(\mathit{x}\right)-\mathcal{M}{\left(g\left(\mathit{x}\right)\right)}_{i_1{i}_2\dots i_n}{\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {\mathit{v}}_{i_n}}}$$

(20)

where n is the number of orthogonal coordinates defined by the vector $\mathit{x}=\{x_1,x_2,\dots ,x_n\}$, ${\mathit{v}}_{i_k}\left(x_k\right)$ is the $i_k$th element of a vector function of the variable $x_k$, $\mathcal{M}$ is an n-dimensional tensor that is a function of the boundary constraints defined in $c\left(\mathit{x}\right)$, and $g\left(\mathit{x}\right)$ is the free-function.

In Equation (20), the term $A\left(\mathit{x}\right)$ represents any function satisfying the boundary constraints defined by $c\left(\mathit{x}\right)$ and the term $B\left(\mathit{x}\right)$ represents all possible functions that are zero on the boundary constraints. The subsections that follow explain how to construct the $\mathcal{M}$ tensor and the ${\mathit{v}}_{i_k}$ vectors for assigned boundary constraints, and provides a proof that the tensor formulation of the ToC defined by Equation (20) satisfies all boundary constraints defined by $c\left(\mathit{x}\right)$, independently of the choice of the free function, $g\left(\mathit{x}\right)$.

Consider a generic boundary constraint on the $x_k=p$ hyperplane, where $k\in [1,n]$. This constraint specifies the d-derivative of the constraint function $c\left(\mathit{x}\right)$ evaluated at $x_k=p$ and it is indicated by ${}^k{c}_p^d:={\left.{\displaystyle \frac{{\partial }^{d}c\left(\mathit{x}\right)}{\partial x_k^d}}\right|}_{x_k=p}$. Consider a set of ${\ell }_k$ constraints defined in various $x_k$ hyperplanes. This set of constraints is indicated by ${}^k{c}_{{\mathit{p}}^k}^{{\mathit{d}}^k}$, where ${\mathit{d}}^k$ and ${\mathit{p}}^k$ are vectors of ${\ell }_k$ elements indicating the order of derivatives and the values of $x_k$ where the boundary constraints are defined, respectively. A specific boundary constraint, e.g. the mth boundary constraint, can then be written as ${}^k{c}_{{\mathit{p}}_m^k}^{{\mathit{d}}_m^k}$.

Additionally, let us define an operator, called the boundary constraint operator, whose purpose is to take the dth derivative with respect to coordinate $x_k$ and then evaluate that function at $x_k=p$. Equation (21) shows the idea.

$${}^k{b}_p^d\left[f\right]\equiv \frac{{\partial }^{d}f}{\partial x_k^d}{|}_{(x_1,\dots ,x_{k-1},p,x_{k+1},\dots ,x_n)}$$

(21)

In general, for a function of n variables, the boundary constraint operator identifies an $n-1$-dimensional manifold. As the boundary constraint operator is used throughout this proof, it is important to note its properties when acting on sums and products of functions. Equation (22) shows how the boundary constraint operator acts on sums, and Equation (23) shows how the boundary constraint operator acts on products.

$${}^k{b}_p^d[f_1+f_2]={}^k{b}_p^d\left[f_1\right]+{}^k{b}_p^d\left[f_2\right]$$

(22)

$${}^k{b}_p^d\left[f_1{f}_2\right]=\left\{\begin{array}{cc}{}^k{b}_p^d\left[f_1\right]{}^k{b}_p^d\left[f_2\right],\hfill & d=0\hfill \\ {}^k{b}_p^d\left[f_1\right]f_2+f_1{}^k{b}_p^d\left[f_2\right],\hfill & d>0\hfill \end{array}\right.$$

(23)

This section shows how to build the $\mathcal{M}$ tensor and the vectors $\mathit{v}$ given the boundary constraints defined by the boundary constraint operators. Moreover, this section contains a proof that, in Equation (20), the boundary constraints defined by $c\left(\mathit{x}\right)$ satisfy the function $A\left(\mathit{x}\right)$ and, by extension, the function $B\left(\mathit{x}\right)$ projects the free-function $g\left(\mathit{x}\right)$ onto the sub-space of functions that are zero on the boundary constraints. Then, it follows that the expression for the ToC surface given in Equation (20) represents all possible functions that meet the boundary defined by the boundary constraint operators.

7.1. The $\mathcal{M}$ Tensor

There is a step-by-step method for constructing the $\mathcal{M}$ tensor.

• The element of $\mathcal{M}$ for all indices equal to 1 is 0 (i.e., ${\mathcal{M}}_{11\dots 1}=0$).
• The first order tensor obtained by keeping the kth dimension’s index and setting all other dimension’s indices to 1 can be written as,

  $${\mathcal{M}}_{1,\dots ,1,i_k,1,\dots ,1}={}^k{c}_{{\mathit{p}}^k}^{{\mathit{d}}^k},\mathrm{where}{i}_k\in [2,{\ell }_k+1],$$

  where the vector ${}^k{c}_{{\mathit{p}}^k}^{{\mathit{d}}^k}$ contains the ${\ell }_k$ boundary constraints specified along the $x_k$-axis. For example, consider the following ${\ell }_7=3$ constraints on the $k=7$th axis,

  $${}^7{c}_{{\mathit{p}}^7}^{{\mathit{d}}^7}:=\left\{{c|}_{x_7=-0.3},{\left.{\displaystyle \frac{{\partial }^{4}c}{\partial x_7^4}}\right|}_{x_7=0.5},{\left.{\displaystyle \frac{\partial c}{\partial x_7}}\right|}_{x_7=1.1}\right\}\mathrm{then}:\left\{\begin{array}{c}{\mathit{d}}^7=\{0,4,1\}\hfill \\ {\mathit{p}}^7=\{-0.3,0.5,1.1\}.\hfill \end{array}\right.$$
• The generic element of the tensor is ${\mathcal{M}}_{i_1{i}_2\dots i_n}$, where at least two indices are different from 1. Let m be the number of indices different from 1. Note that m is also the number of constraint “intersections”. In this case, the generic element of the $\mathcal{M}$ tensor is provided by,

  $${\mathcal{M}}_{i_1{i}_2\dots i_n}={}^1{b}_{{\mathit{p}}_{i_1-1}^1}^{{\mathit{d}}_{i_1-1}^1}\left[{}^2{b}_{{\mathit{p}}_{i_2-1}^2}^{{\mathit{d}}_{i_2-1}^2}\left[\dots \left[{}^n{b}_{{\mathit{p}}_{i_n-1}^n}^{{\mathit{d}}_{i_n-1}^n}\left[c\left(\mathit{x}\right)\right]\right]\dots \right]\right]{(-1)}^{m+1}.$$

  (24)
  If $c\left(\mathit{x}\right)\in C^s$, where ${\displaystyle s=\sum _{k=1}^n{\mathit{d}}_{i_k-1}^k}$, then Clairaut’s theorem states that the sequence of boundary constraint operators provided in Equation (24) can be freely permutated. This permutation becomes obvious by multiple applications of the theorem. For example,

  $$f_{xyy}={\left(f_{xy}\right)}_y={\left(f_{yx}\right)}_y={\left(f_y\right)}_{xy}={\left(f_y\right)}_{yx}=f_{yyx}.$$

To better clarify how to use Equation (24), consider the example of the following constraints in three-dimensional space.

$${c\left(\mathit{x}\right)|}_{x_1=0}{,c\left(\mathit{x}\right)|}_{x_1=1},c\left(\mathit{x}\right){{|}_{x_2=0},{\displaystyle \frac{\partial c\left(\mathit{x}\right)}{\partial x_2}}{|}_{x_2=0},c\left(\mathit{x}\right)|}_{x_3=0},\mathrm{and}{\displaystyle \frac{\partial c\left(\mathit{x}\right)}{\partial x_3}}{|}_{x_3=0}$$

• From Step 1: $M_{111}=0$
• From Step 2:

  $$\begin{array}{cc}\hfill M_{i_{1}11}& =\left\{\begin{array}{ccc}0,& c(0,x_2,x_3),& c(1,x_2,x_3)\end{array}\right\}=\left\{\begin{array}{ccc}0,& {}^1{b}_0^0\left[c\left(\mathit{x}\right)\right],& {}^1{b}_1^0\left[c\left(\mathit{x}\right)\right]\end{array}\right\}\hfill \\ \hfill M_{1i_{2}1}& =\left\{\begin{array}{ccc}0,& c(x_1,0,x_3),& {\displaystyle \frac{\partial c}{\partial x_2}}(x_1,0,x_3)\end{array}\right\}=\left\{\begin{array}{ccc}0,& {}^2{b}_0^0\left[c\left(\mathit{x}\right)\right],& {}^2{b}_0^1\left[c\left(\mathit{x}\right)\right]\end{array}\right\}\hfill \\ \hfill M_{11i_3}& =\left\{\begin{array}{ccc}0,& c(x_1,x_3,0),& {\displaystyle \frac{\partial c}{\partial x_3}}(x_1,x_2,0)\end{array}\right\}=\left\{\begin{array}{ccc}0,& {}^3{b}_0^0\left[c\left(\mathit{x}\right)\right],& {}^3{b}_0^1\left[c\left(\mathit{x}\right)\right]\end{array}\right\}\hfill \end{array}$$
• From Step 3, a single example is provided,

  $${\mathcal{M}}_{323}={}^1{b}_1^0\left[{}^2{b}_0^0\left[{}^3{c}_0^1\left(\mathit{x}\right)\right]\right]{(-1)}^4={\displaystyle \frac{\partial c\left(\mathit{x}\right)}{\partial x_3}}{|}_{\begin{array}{c}{x}_1=1\\ x_2=0\\ x_3=0\end{array}}$$

  which, thanks to Clairaut’s theorem, can also be written as,

  $${\mathcal{M}}_{323}={}^2{b}_0^0\left[{}^3{b}_0^1\left[{}^1{c}_1^0\right]\right]{(-1)}^4={}^3{b}_0^1\left[{}^1{b}_1^0\left[{}^2{c}_0^0\right]\right]{(-1)}^4.$$
  Three additional examples are given to help further illustrate the procedure,

  $$M_{132}=-{\left.{\displaystyle \frac{\partial c\left(\mathit{x}\right)}{\partial x_2}}\right|}_{\begin{array}{c}{x}_2=0\\ x_3=0\end{array}},M_{221}=-c(0,0,x_3),\mathrm{and}{M}_{333}={\displaystyle \frac{{\partial }^{2}c\left(\mathit{x}\right)}{\partial x_2\partial x_3}}{|}_{\begin{array}{c}{x}_1=1\\ x_2=0\\ x_3=0\end{array}}$$

7.2. The v Vectors

Each vector, ${\mathit{v}}_k$, is associated with the ${\ell }_k$ constraints that are specified by ${}^k{c}_{{\mathit{p}}^k}^{{\mathit{d}}^k}$. The ${\mathit{v}}_k$ vector is built as follows,

$${\mathit{v}}_k=\left\{\begin{array}{ccccc}1,& {\displaystyle \sum _{i=1}^{{\ell }_k}{\alpha }_{i1}{h}_i\left(x_k\right),}& {\displaystyle \sum _{i=1}^{{\ell }_k}{\alpha }_{i2}{h}_i\left(x_k\right),}& \dots ,& {\displaystyle \sum _{i=1}^{{\ell }_k}{\alpha }_{i{\ell }_k}{h}_i\left(x_k\right)}\end{array}\right\},$$

where $h_i\left(x_k\right)$ are ${\ell }_k$ linearly independent functions. The simplest set of linearly independent functions are monomials, that is, $h_i\left(x_k\right)=x_k^{i-1}$. The ${\ell }_k\times {\ell }_k$ coefficients, ${\alpha }_{ij}$, can be computed by matrix inversion,

$$\left[\begin{array}{cccc}{}^k{b}_{p_1}^{d_1}\left[h_1\right]& {}^k{b}_{p_1}^{d_1}\left[h_2\right]& \dots & {}^k{b}_{p_1}^{d_1}\left[h_{{\ell }_k}\right]\\ {}^k{b}_{p_2}^{d_2}\left[h_1\right]& {}^k{b}_{p_2}^{d_2}\left[h_2\right]& \dots & {}^k{b}_{p_2}^{d_2}\left[h_{{\ell }_k}\right]\\ ⋮& ⋮& \ddots & ⋮\\ {}^k{b}_{p_{{\ell }_k}}^{d_{{\ell }_k}}\left[h_1\right]& {}^k{b}_{p_{{\ell }_k}}^{d_{{\ell }_k}}\left[h_2\right]& \dots & {}^k{b}_{p_{{\ell }_k}}^{d_{{\ell }_k}}\left[h_{{\ell }_k}\right]\end{array}\right]\left[\begin{array}{cccc}{\alpha }_{11}& {\alpha }_{12}& \dots & {\alpha }_{1{\ell }_k}\\ {\alpha }_{21}& {\alpha }_{22}& \dots & {\alpha }_{2{\ell }_k}\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{{\ell }_{k}1}& {\alpha }_{{\ell }_{k}2}& \dots & {\alpha }_{{\ell }_k{\ell }_k}\end{array}\right]=\left[\begin{array}{cccc}1& 0& \dots & 0\\ 0& 1& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1\end{array}\right].$$

(25)

To supplement the above explanation, let us look at the example of Dirichlet boundary conditions on $x_1$ from the example in Section 7.1. There are two boundary conditions, ${c\left(\mathit{x}\right)|}_{x_1=0}$ and ${c\left(\mathit{x}\right)|}_{x_1=1}$, and thus two linearly independent functions are needed,

$$\begin{array}{c}\hfill {\mathit{v}}_{i_1}=\left\{\begin{array}{ccc}1,& {\alpha }_{11}{h}_1\left(x_1\right)+{\alpha }_{21}{h}_2\left(x_1\right),& {\alpha }_{12}{h}_1\left(x_1\right)+{\alpha }_{22}{h}_2\left(x_1\right)\end{array}\right\}.\end{array}$$

Let us consider, $h_1\left(x_1\right)=1$ and $h_2\left(x_1\right)=x_1$. Then, following Equation (25),

$$\left[\begin{array}{cc}{}^1{b}_0^0\left[1\right]& {}^1{b}_0^0\left[x\right]\\ {}^2{b}_1^0\left[1\right]& {}^2{b}_1^0\left[x\right]\end{array}\right]\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\to \left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\begin{array}{cc}\hfill \left[r\right]1& 0\\ \hfill -1& 1\end{array},$$

and substituting the values of ${\alpha }_{ij}$, we obtain ${\mathit{v}}_{i_1}=\left\{\begin{array}{ccc}1,& 1-x_1,& x_1\end{array}\right\}$.

7.3. Proof

This section demonstrates that the term $A\left(\mathit{x}\right)$ from Equation (20) generates a surface satisfying the boundary constraints defined by the function $c\left(\mathit{x}\right)$. First, it is shown that $A\left(\mathit{x}\right)$ satisfies boundary constraints on the value, and then that $A\left(\mathit{x}\right)$ satisfies boundary constraints on arbitrary-order derivatives.

Equation (23) for $d=0$ allows us to write,

$${}^k{b}_{p_{q-1}}^0\left[A\left(\mathit{x}\right)\right]={}^k{b}_{p_{q-1}}^0\left[{\mathcal{M}}_{i_1{i}_2\dots i_k\dots i_n}\right]{\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {}^k{b}_{p_{q-1}}^0\left[{\mathit{v}}_{i_k}\right]\dots {\mathit{v}}_{i_n}.$$

(26)

The boundary constraint operator applied to ${\mathit{v}}_k$ yields,

$${}^k{b}_{p_{q-1}}^0\left[{\mathit{v}}_{i_k}\right]=\left\{\begin{array}{cc}=1,\hfill & i_k=1,q\hfill \\ =0,\hfill & i_k\ne 1,q.\hfill \end{array}\right.$$

(27)

Since the only nonzero terms are associated with $i_k=1,q$, we have,

$${}^k{b}_{p_{q-1}}^0\left[A\left(\mathit{x}\right)\right]=\left({}^k{b}_{p_{q-1}}^0\left[{\mathcal{M}}_{i_1{i}_2\dots 1\dots i_n}\right]+{}^k{b}_{p_{q-1}}^0\left[{\mathcal{M}}_{i_1{i}_2\dots q\dots i_n}\right]\right){\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {\mathit{v}}_{i_n}.$$

(28)

Applying the boundary constraint operator to the $n-1$-dimensional $\mathcal{M}$ tensor where index $i_k=q$ has no effect, because all of the functions already have coordinate $x_k$ substituted for the value $p_{q-1}$ (see Equation (24)). Moreover, applying the boundary constraint operator to the $\mathcal{M}$ tensor where index $i_k=1$ causes all terms in the sum within the parenthesis in Equation (28) to cancel each other, except when all of the non-$i_k$ indices are equal to one. This leads to Equation (29).

$${}^k{b}_{p_{q-1}}^0\left[A\left(\mathit{x}\right)\right]=\left({\mathcal{M}}_{11\dots 1\dots 1}+{\mathcal{M}}_{11\dots q\dots 1}\right){\mathit{v}}_1{\mathit{v}}_1\dots {\mathit{v}}_1$$

(29)

Since ${\mathit{v}}_j=1$ when $j=1$ and ${\mathcal{M}}_{11\dots 1}=0$ by definition, then,

$${}^k{b}_{p_{q-1}}^0\left[A\left(\mathit{x}\right)\right]={\mathcal{M}}_{11\dots q\dots 1}=c(x_1,x_2,\dots ,p_{q-1},\dots ,x_n),$$

which proves Equation (20) works for boundary constraints on the value.

Now, we show that Equation (20) holds for arbitrary-order derivative type boundary constraints. Equation (23) for $d>0$ allows us to write,

$${}^k{b}_{p_{q-1}}^{d_{q-1}}\left[A\left(\mathit{x}\right)\right]={}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathcal{M}}_{i_1{i}_2\dots i_k\dots i_n}\right]{\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {\mathit{v}}_{i_k}\dots {\mathit{v}}_{i_n}+{\mathcal{M}}_{i_1{i}_2\dots i_k\dots i_n}{\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathit{v}}_{i_k}\right]\dots {\mathit{v}}_{i_n}.$$

(30)

From Equation (23), we note that boundary constraint operators that take a derivative follow the usual product rule when applied to a product. Moreover, we note that all of the $\mathit{v}$ vectors except ${\mathit{v}}_{i_k}$ do not depend on $x_k$, thus applying the boundary constraint operator to them results in a vector of zeros. Applying the boundary constraint operator to ${\mathit{v}}_{i_k}$ yields,

$${}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathit{v}}_{i_k}\right]=\left\{\begin{array}{cc}=1,\hfill & i_k=q\hfill \\ =0,\hfill & i_k\ne q,\hfill \end{array}\right.$$

and applying the boundary constraint operator to $\mathcal{M}$ yields,

$${}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathcal{M}}_{i_1{i}_2\dots 1\dots i_n}\right]=\left\{\begin{array}{cc}={}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathcal{M}}_{i_1{i}_2\dots 1\dots i_n}\right],\hfill & i_k=1\hfill \\ =0,\hfill & i_k\ne 1.\hfill \end{array}\right.$$

Substituting these simplifications into $A\left(\mathit{x}\right)={\mathcal{M}}_{i_1{i}_2\dots i_k\dots i_n}{\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {\mathit{v}}_{i_k}\dots {\mathit{v}}_{i_n}$, after applying the boundary constraint operator, results in Equation (31).

$${}^k{b}_{p_{q-1}}^{d_{q-1}}\left[A\left(\mathit{x}\right)\right]=\left({}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathcal{M}}_{i_1{i}_2\dots 1\dots i_n}\right]+{\mathcal{M}}_{i_1{i}_2\dots q\dots i_n}\right){\mathit{v}}_{i_1}{\mathit{v}}_{i_2}\dots {\mathit{v}}_{i_n}$$

(31)

Similar to the proof for value-based boundary constraints, based on Equation (24), all terms in the sum within the parenthesis in Equation (31) cancel each other, except when all of the non-$i_k$ indices are equal to one. Thus, Equation (31) can be simplified to Equation (32).

$${}^k{b}_{p_{q-1}}^{d_{q-1}}\left[A\left(\mathit{x}\right)\right]=\left({}^k{b}_{p_{q-1}}^{d_{q-1}}\left[{\mathcal{M}}_{11\dots 1\dots 1}\right]+{\mathcal{M}}_{11\dots q\dots 1}\right){\mathit{v}}_1{\mathit{v}}_1\dots {\mathit{v}}_1$$

(32)

Again, all of the vectors $\mathit{v}$ were designed such that their first component is 1, and the value of the element of $\mathcal{M}$ for all indices equal to 1 is 0. Therefore, Equation (32) simplifies to,

$${}^k{b}_{p_{q-1}}^{d_{q-1}}\left[A\left(\mathit{x}\right)\right]={\mathcal{M}}_{11\dots q\dots 1}={\displaystyle \frac{{\partial }^{d}c\left(\mathit{x}\right)}{\partial x_k^d}}{|}_{x_k=p_{q-1}},$$

which proves Equation (20) works for arbitrary-order derivative boundary constraints.

In conclusion, the term $A\left(\mathit{x}\right)$ from Equation (20) generates a manifold satisfying the boundary constraints given in terms of arbitrary-order derivative in n-dimensional space. The term $B\left(\mathit{x}\right)$ from Equation (20) projects any free function $g\left(\mathit{x}\right)$ onto the space of functions that are vanishing at the specified boundary constraints. As a result, Equation (20) can be used to produce the family of all possible functions satisfying assigned boundary constraints (functions or derivatives) in rectangular domains in n-dimensional space.

8. Conclusions

This paper extends to n-dimensional spaces the Univariate Theory of Connections (ToC), introduced in Ref. [1]. First, it provides a mathematical tool to express all possible surfaces subject to constraint functions and arbitrary-order derivatives in a boundary rectangular domain, and then it extends the results to the multivariate case by providing the Multivariate Theory of Connections, which allows one to obtain n-dimensional manifolds subject to any-order derivative boundary constraints.

In particular, if the constraints are provided along one axis only, then this paper shows that the univariate ToC, as defined in Ref. [1], can be adopted to describe all possible surfaces satisfying the constraints. If the boundary constraints are defined in a rectangular domain, then the constrained expression is found in the form $f\left(\mathit{x}\right)=A\left(\mathit{x}\right)+B\left(\mathit{x}\right)$, where $A\left(\mathit{x}\right)$ can be any function satisfying the constraints and $B\left(\mathit{x}\right)$ describes all functions that are vanishing at the constraints. This is obtained by introducing a free function, $g\left(\mathit{x}\right)$, into the function $B\left(\mathit{x}\right)$ in such a way that $B\left(\mathit{x}\right)$ is zero at the constraints no matter what the $g\left(\mathit{x}\right)$ is. This way, by spanning all possible $g\left(\mathit{x}\right)$ surfaces (even discontinuous, null, or piece-wise defined) the resulting $B\left(\mathit{x}\right)$ generates all surfaces that are zero at the constraints and, consequently, $f\left(\mathit{x}\right)=A\left(\mathit{x}\right)+B\left(\mathit{x}\right)$, describes all surfaces satisfying the constraints defined in the rectangular boundary domain. The function $A\left(\mathit{x}\right)$ has been selected as a Coons surface [11] and, in particular, a Coons surface is obtained if $g\left(\mathit{x}\right)=0$ is selected. All possible combinations of Dirichlet and Neumann constraints are also provided in Appendix A.

The last section provides the Multivariate Theory of Connections extension, which is a mathematical tool to transform n-dimensional constraint optimization problems subject to constraints on the boundary value and any-order derivative into unconstrained optimization problems. The number of applications of the Multivariate Theory of Connections are many, especially in the area of partial and stochastic differential equations: the main subjects of our current research.