All the concepts and algorithms developed for approximation of 1D functions \(f(x)\) can readily be extended to 2D functions \(f(x,y)\) and 3D functions \(f(x,y,z)\). Basically, the extensions consists of defining basis functions \({\psi}_i(x,y)\) or \({\psi}_i(x,y,z)\) over some domain \(\Omega\), and for the least squares and Galerkin methods, the integration is done over \(\Omega\).

As in 1D, the least squares and projection/Galerkin methods two lead to linear systems

\[\begin{split}\sum_{j\in{\mathcal{I}_s}} A_{i,j}c_j &= b_i,\quad i\in{\mathcal{I}_s},\\
A_{i,j} &= ({\psi}_i,{\psi}_j),\\
b_i &= (f,{\psi}_i),\end{split}\]

where the inner product of two functions \(f(x,y)\) and \(g(x,y)\) is defined
completely analogously to the 1D case *(4.2)*:

\[(f,g) = \int_\Omega f(x,y)g(x,y) dx dy\]

One straightforward way to construct a basis in 2D is to combine 1D basis functions. Say we have the 1D vector space

(1)\[ V_x = \mbox{span}\{ \hat{\psi}_0(x),\ldots,\hat{\psi}_{N_x}(x)\}\]\[ {\thinspace .}\]

A similar space for variation in \(y\) can be defined,

(2)\[ V_y = \mbox{span}\{ \hat{\psi}_0(y),\ldots,\hat{\psi}_{N_y}(y)\}\]\[ {\thinspace .}\]

We can then form 2D basis functions as *tensor products* of 1D basis functions.

Tensor products

Given two vectors \(a=(a_0,\ldots,a_M)\) and \(b=(b_0,\ldots,b_N)\),
their *outer tensor product*, also called the *dyadic product*,
is \(p=a\otimes b\), defined through

\[p_{i,j}=a_ib_j,\quad i=0,\ldots,M,\ j=0,\ldots,N{\thinspace .}\]

In the tensor terminology,
\(a\) and \(b\) are first-order tensors (vectors with one index, also termed
rank-1 tensors), and then their outer
tensor product is a second-order tensor (matrix with two indices, also
termed rank-2 tensor). The
corresponding *inner tensor product* is the well-known scalar or dot
product of two vectors: \(p=a\cdot b = \sum_{j=0}^N a_jb_j\). Now,
\(p\) is a rank-0 tensor.

Tensors are typically represented by arrays in computer code. In the above example, \(a\) and \(b\) are represented by one-dimensional arrays of length \(M\) and \(N\), respectively, while \(p=a\otimes b\) must be represented by a two-dimensional array of size \(M\times N\).

Tensor products can be used in a variety of context.

Given the vector spaces \(V_x\) and \(V_y\) as defined in (1) and (2), the tensor product space \(V=V_x\otimes V_y\) has a basis formed as the tensor product of the basis for \(V_x\) and \(V_y\). That is, if \(\left\{ {\varphi}_i(x) \right\}_{i\in{\mathcal{I}_x}}\) and \(\left\{ {\varphi}_i(y) \right\}_{i\in {\mathcal{I}_y}}\) are basis for \(V_x\) and \(V_y\), respectively, the elements in the basis for \(V\) arise from the tensor product: \(\left\{ {\varphi}_i(x){\varphi}_j(y) \right\}_{i\in {\mathcal{I}_x},j\in {\mathcal{I}_y}}\). The index sets are \(I_x=\{0,\ldots,N_x\}\) and \(I_y=\{0,\ldots,N_y\}\).

The notation for a basis function in 2D can employ a double index as in

\[{\psi}_{p,q}(x,y) = \hat{\psi}_p(x)\hat{\psi}_q(y),
\quad p\in{\mathcal{I}_x},q\in{\mathcal{I}_y}{\thinspace .}\]

The expansion for \(u\) is then written as a double sum

\[u = \sum_{p\in{\mathcal{I}_x}}\sum_{q\in{\mathcal{I}_y}} c_{p,q}{\psi}_{p,q}(x,y){\thinspace .}\]

Alternatively, we may employ a single index,

\[{\psi}_i(x,y) = \hat{\psi}_p(x)\hat{\psi}_q(y),\]

and use the standard form for \(u\),

\[u = \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x,y){\thinspace .}\]

The single index is related to the double index through \(i=p (N_y+1) + q\) or \(i=q (N_x+1) + p\).

Suppose we choose \(\hat{\psi}_p(x)=x^p\), and try an approximation with \(N_x=N_y=1\):

\[{\psi}_{0,0}=1,\quad {\psi}_{1,0}=x, \quad {\psi}_{0,1}=y,
\quad {\psi}_{1,1}=xy
{\thinspace .}\]

Using a mapping to one index like \(i=q (N_x+1) + p\), we get

\[{\psi}_0=1,\quad {\psi}_1=x, \quad {\psi}_2=y,\quad{\psi}_3 =xy
{\thinspace .}\]

With the specific choice \(f(x,y) = (1+x^2)(1+2y^2)\) on \(\Omega = [0,L_x]\times [0,L_y]\), we can perform actual calculations:

\[\begin{split}A_{0,0} &= ({\psi}_0,{\psi}_0) = \int_0^{L_y}\int_{0}^{L_x}
{\psi}_0(x,y)^2 dx dy = \int_0^{L_y}\int_{0}^{L_x}dx dy = L_xL_y,\\
A_{1,0} &= ({\psi}_1,{\psi}_0) = \int_0^{L_y}\int_{0}^{L_x} x dxdy =
{\frac{1}{2}}L_x^2L_y,\\
A_{0,1} &= ({\psi}_0,{\psi}_1) = \int_0^{L_y}\int_{0}^{L_x} y dxdy =
{\frac{1}{2}}L_y^2L_x,\\
A_{0,1} &= ({\psi}_0,{\psi}_1) = \int_0^{L_y}\int_{0}^{L_x} xy dxdy =
\int_0^{L_y}ydy \int_{0}^{L_x} xdx =
\frac{1}{4}L_y^2L_x^2
{\thinspace .}\end{split}\]

The right-hand side vector has the entries

\[\begin{split}b_{0} &= ({\psi}_0,f) = \int_0^{L_y}\int_{0}^{L_x}1\cdot (1+x^2)(1+2y^2) dxdy\\
&= \int_0^{L_y}(1+2y^2)dy \int_{0}^{L_x} (1+x^2)dx
= (L_y + \frac{2}{3}L_y^3)(L_x + \frac{1}{3}L_x^3)\\
b_{1} &= ({\psi}_1,f) = \int_0^{L_y}\int_{0}^{L_x} x(1+x^2)(1+2y^2) dxdy\\
&=\int_0^{L_y}(1+2y^2)dy \int_{0}^{L_x} x(1+x^2)dx
= (L_y + \frac{2}{3}L_y^3)({\frac{1}{2}}L_x^2 + \frac{1}{4}L_x^4)\\
b_{2} &= ({\psi}_2,f) = \int_0^{L_y}\int_{0}^{L_x} y(1+x^2)(1+2y^2) dxdy\\
&= \int_0^{L_y}y(1+2y^2)dy \int_{0}^{L_x} (1+x^2)dx
= ({\frac{1}{2}}L_y + {\frac{1}{2}}L_y^4)(L_x + \frac{1}{3}L_x^3)\\
b_{3} &= ({\psi}_2,f) = \int_0^{L_y}\int_{0}^{L_x} xy(1+x^2)(1+2y^2) dxdy\\
&= \int_0^{L_y}y(1+2y^2)dy \int_{0}^{L_x} x(1+x^2)dx
= ({\frac{1}{2}}L_y^2 + {\frac{1}{2}}L_y^4)({\frac{1}{2}}L_x^2 + \frac{1}{4}L_x^4)
{\thinspace .}\end{split}\]

There is a general pattern in these calculations that we can explore. An arbitrary matrix entry has the formula

\[\begin{split}A_{i,j} &= ({\psi}_i,{\psi}_j) = \int_0^{L_y}\int_{0}^{L_x}
{\psi}_i{\psi}_j dx dy \\
&= \int_0^{L_y}\int_{0}^{L_x}
{\psi}_{p,q}{\psi}_{r,s} dx dy
= \int_0^{L_y}\int_{0}^{L_x}
\hat{\psi}_p(x)\hat{\psi}_q(y)\hat{\psi}_r(x)\hat{\psi}_s(y) dx dy\\
&= \int_0^{L_y} \hat{\psi}_q(y)\hat{\psi}_s(y)dy
\int_{0}^{L_x} \hat{\psi}_p(x) \hat{\psi}_r(x) dx\\
&= \hat A^{(x)}_{p,r}\hat A^{(y)}_{q,s},\end{split}\]

where

\[\hat A^{(x)}_{p,r} = \int_{0}^{L_x} \hat{\psi}_p(x) \hat{\psi}_r(x) dx,
\quad
\hat A^{(y)}_{q,s} = \int_0^{L_y} \hat{\psi}_q(y)\hat{\psi}_s(y)dy,\]

are matrix entries for one-dimensional approximations. Moreover, \(i=q N_y+q\) and \(j=s N_y+r\).

With \(\hat{\psi}_p(x)=x^p\) we have

\[\hat A^{(x)}_{p,r} = \frac{1}{p+r+1}L_x^{p+r+1},\quad
\hat A^{(y)}_{q,s} = \frac{1}{q+s+1}L_y^{q+s+1},\]

and

\[A_{i,j} = \hat A^{(x)}_{p,r} \hat A^{(y)}_{q,s} =
\frac{1}{p+r+1}L_x^{p+r+1} \frac{1}{q+s+1}L_y^{q+s+1},\]

for \(p,r\in{\mathcal{I}_x}\) and \(q,s\in{\mathcal{I}_y}\).

Corresponding reasoning for the right-hand side leads to

\[\begin{split}b_i &= ({\psi}_i,f) = \int_0^{L_y}\int_{0}^{L_x}{\psi}_i f\,dxdx\\
&= \int_0^{L_y}\int_{0}^{L_x}\hat{\psi}_p(x)\hat{\psi}_q(y) f\,dxdx\\
&= \int_0^{L_y}\hat{\psi}_q(y) (1+2y^2)dy
\int_0^{L_y}\hat{\psi}_p(x) x^p (1+x^2)dx\\
&= \int_0^{L_y} y^q (1+2y^2)dy
\int_0^{L_y}x^p (1+x^2)dx\\
&= (\frac{1}{q+1} L_y^{q+1} + \frac{2}{q+3}L_y^{q+3})
(\frac{1}{p+1} L_x^{p+1} + \frac{2}{q+3}L_x^{p+3})\end{split}\]

Choosing \(L_x=L_y=2\), we have

\[\begin{split}A =
\left[\begin{array}{cccc}
4 & 4 & 4 & 4\\
4 & \frac{16}{3} & 4 & \frac{16}{3}\\
4 & 4 & \frac{16}{3} & \frac{16}{3}\\
4 & \frac{16}{3} & \frac{16}{3} & \frac{64}{9}
\end{array}\right],\quad
b = \left[\begin{array}{c}
\frac{308}{9}\\\frac{140}{3}\\44\\60\end{array}\right],
\quad c = \left[
\begin{array}{r}
-\frac{1}{9} \\
\frac{4}{3} \\
- \frac{2}{3} \\
8
\end{array}\right]
{\thinspace .}\end{split}\]

Figure *Approximation of a 2D quadratic function (left) by a 2D bilinear function (right) using the Galerkin or least squares method* illustrates the result.

The `least_squares` function from
the section *Orthogonal basis functions* and/or the
file approx1D.py
can with very small modifications solve 2D approximation problems.
First, let `Omega` now be a list of the intervals in \(x\) and \(y\) direction.
For example, \(\Omega = [0,L_x]\times [0,L_y]\) can be represented
by `Omega = [[0, L_x], [0, L_y]]`.

Second, the symbolic integration must be extended to 2D:

```
import sympy as sp
integrand = psi[i]*psi[j]
I = sp.integrate(integrand,
(x, Omega[0][0], Omega[0][1]),
(y, Omega[1][0], Omega[1][1]))
```

provided `integrand` is an expression involving the `sympy` symbols `x`
and `y`.
The 2D version of numerical integration becomes

```
if isinstance(I, sp.Integral):
integrand = sp.lambdify([x,y], integrand)
I = sp.mpmath.quad(integrand,
[Omega[0][0], Omega[0][1]],
[Omega[1][0], Omega[1][1]])
```

The right-hand side integrals are modified in a similar way.

Third, we must construct a list of 2D basis functions. Here are two examples based on tensor products of 1D “Taylor-style” polynomials \(x^i\) and 1D sine functions \(\sin((i+1)\pi x)\):

```
def taylor(x, y, Nx, Ny):
return [x**i*y**j for i in range(Nx+1) for j in range(Ny+1)]
def sines(x, y, Nx, Ny):
return [sp.sin(sp.pi*(i+1)*x)*sp.sin(sp.pi*(j+1)*y)
for i in range(Nx+1) for j in range(Ny+1)]
```

The complete code appears in approx2D.py.

The previous hand calculation where a quadratic \(f\) was approximated by a bilinear function can be computed symbolically by

```
>>> from approx2D import *
>>> f = (1+x**2)*(1+2*y**2)
>>> psi = taylor(x, y, 1, 1)
>>> Omega = [[0, 2], [0, 2]]
>>> u, c = least_squares(f, psi, Omega)
>>> print u
8*x*y - 2*x/3 + 4*y/3 - 1/9
>>> print sp.expand(f)
2*x**2*y**2 + x**2 + 2*y**2 + 1
```

We may continue with adding higher powers to the basis:

```
>>> psi = taylor(x, y, 2, 2)
>>> u, c = least_squares(f, psi, Omega)
>>> print u
2*x**2*y**2 + x**2 + 2*y**2 + 1
>>> print u-f
0
```

For \(N_x\geq 2\) and \(N_y\geq 2\) we recover the exact function \(f\), as
expected, since in that case \(f\in V\) (see
the section *Perfect approximation*).

Extension to 3D is in principle straightforward once the 2D extension is understood. The only major difference is that we need the repeated outer tensor product,

\[V = V_x\otimes V_y\otimes V_z{\thinspace .}\]

In general, given vectors (first-order tensors) \(a^{(q)} = (a^{(q)}_0,\ldots,a^{(q)}_{N_q}\), \(q=0,\ldots,m\), the tensor product \(p=a^{(0)}\otimes\cdots\otimes a^{m}\) has elements

\[p_{i_0,i_1,\ldots,i_m} = a^{(0)}_{i_1}a^{(1)}_{i_1}\cdots a^{(m)}_{i_m}{\thinspace .}\]

The basis functions in 3D are then

\[{\psi}_{p,q,r}(x,y,z) = \hat{\psi}_p(x)\hat{\psi}_q(y)\hat{\psi}_r(z),\]

with \(p\in{\mathcal{I}_x}\), \(q\in{\mathcal{I}_y}\), \(r\in{\mathcal{I}_z}\). The expansion of \(u\) becomes

\[u(x,y,z) = \sum_{p\in{\mathcal{I}_x}}\sum_{q\in{\mathcal{I}_y}}\sum_{r\in{\mathcal{I}_z}} c_{p,q,r}
{\psi}_{p,q,r}(x,y,z){\thinspace .}\]

A single index can be introduced also here, e.g., \(i=N_xN_yr + q_Nx + p\), \(u=\sum_i c_i{\psi}_i(x,y,z)\).

Use of tensor product spaces

Constructing a multi-dimensional space and basis from tensor products of 1D spaces is a standard technique when working with global basis functions. In the world of finite elements, constructing basis functions by tensor products is much used on quadrilateral and hexahedra cell shapes, but not on triangles and tetrahedra. Also, the global finite element basis functions are almost exclusively denoted by a single index and not by the natural tuple of indices that arises from tensor products.