13.1. Tensor-product discretizations

As you learned when starting double integration in vector calculus, the simplest extension of an interval to two dimensions is a rectangle. We will use a particular notation for rectangles:

(13.1.1)\[ [a,b] \times [c,d] = \bigl\{ (x,y)\in\mathbb{R}^2 : a\le x \le b,\; c\le y \le d \bigr\}.\]

The \(\times\) in this notation is called a tensor product, and a rectangle is the fundamental example of a tensor-product domain. The implication of the tensor product is that each variable independently varies over a fixed set. The simplest three-dimensional tensor-product domain is the cuboid \([a,b]\times[c,d]\times[e,f]\). When the interval is the same in each dimension (that is, the region is a square or a cube), we may write \([a,b]^2\) or \([a,b]^3\). We will limit our discussion to two dimensions henceforth.

The discretization of a two-dimensional tensor-product domain is straightforward.

Definition 13.1.1 :  Tensor-product grid

Given discretizations of two intervals,

(13.1.2)\[a= x_0< x_1 < \cdots < x_m = b, \qquad c = y_0 < y_1 < \cdots < y_n = d, \]

then a tensor-product grid on \([a,b]\times[c,d]\) is the set

(13.1.3)\[ \bigl\{ (x_i,y_j): i=0,\ldots,m,\; j=0,\ldots,n \bigr\}.\]

Functions on grids

The double indexing of the grid set (13.1.3) implies an irresistible connection to matrices. Corresponding to any function \(f(x,y)\) defined on the rectangle is an \((m+1)\times(n+1)\) matrix \(\mathbf{F}\) defined by collecting the values of \(f\) at the points in the grid. This transformation of a function to a matrix is so important that we give it a formal name:

(13.1.4)\[\begin{split}\mathbf{F} = \mtx(f) = \Bigl[f(x_i,y_j)\Bigr]_{\substack{i=0,\ldots,m\\j=0,\ldots,n}}.\end{split}\]


There is potential for confusion because the first dimension of a matrix varies in the vertical direction, while the first coordinate \(x\) varies horizontally. In fact, the Julia plotting routines we use expect the transpose of this arrangement, so that \(x\) varies along columns and \(y\) along rows.

Example 13.1.2

Let the interval \([0,2]\) be divided into \(m=4\) equally sized pieces, and let \([1,3]\) be discretized in \(n=2\) equal pieces. Then the grid in the rectangle \([0,2]\times[1,3]\) is given by all points \((i/2,1+j)\) for all choices \(i=0,1,2,3,4\) and \(j=0,1,2\). If \(f(x,y)=\sin(\pi xy)\), then

\[\begin{split} \mtx(f) = \begin{bmatrix} \sin(\pi\cdot 0\cdot 1) & \sin(\pi\cdot0\cdot 2) & \sin(\pi\cdot0\cdot 3) \\[1mm] \sin\left(\pi\cdot\tfrac{1}{2} \cdot 1 \right) & \sin\left(\pi\cdot\tfrac{1}{2} \cdot 2 \right) & \sin\left(\pi\cdot\tfrac{1}{2} \cdot 3 \right) \\[1mm] \sin\left(\pi \cdot 1 \cdot 1 \right) & \sin\left(\pi \cdot 1 \cdot 2 \right) & \sin\left(\pi \cdot 1 \cdot 3 \right) \\[1mm] \sin\left(\pi\cdot \tfrac{3}{2} \cdot 1 \right) & \sin\left(\pi\cdot\tfrac{3}{2} \cdot 2 \right) & \sin\left(\pi\cdot\tfrac{3}{2} \cdot 3 \right) \\[1mm] \sin\left(\pi \cdot 2 \cdot 1 \right) & \sin\left(\pi \cdot 2 \cdot 2 \right) & \sin\left(\pi \cdot 2 \cdot 3 \right) \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & -1 \\ 0 & 0 & 0 \\ -1 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}. \end{split}\]
Demo 13.1.3

Here is the grid from Example 13.1.2.

m = 4;   x = range(0,2,length=m+1);
n = 2;   y = range(1,3,length=n+1);

For a given \(f(x,y)\) we can find \(\operatorname{mtx}(f)\) by using a comprehension syntax.

f = (x,y) -> cos(π*x*y-y)
F = [ f(x,y) for x in x, y in y ]
5×3 Matrix{Float64}:
  0.540302  -0.416147  -0.989992
  0.841471   0.416147  -0.14112
 -0.540302  -0.416147   0.989992
 -0.841471   0.416147   0.14112
  0.540302  -0.416147  -0.989992

The plots of this section look better using a different graphics engine on the back end:

plotlyjs();  # use better 3D renderer
ArgumentError: Package PlotlyJS not found in current path:
- Run `import Pkg; Pkg.add("PlotlyJS")` to install the PlotlyJS package.

  [1] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:967
  [2] top-level scope
    @ ~/.julia/packages/Plots/8K4be/src/backends.jl:663
  [3] eval
    @ ./boot.jl:373 [inlined]
  [4] _initialize_backend(pkg::Plots.PlotlyJSBackend)
    @ Plots ~/.julia/packages/Plots/8K4be/src/backends.jl:662
  [5] backend(pkg::Plots.PlotlyJSBackend)
    @ Plots ~/.julia/packages/Plots/8K4be/src/backends.jl:176
  [6] #plotlyjs#248
    @ ~/.julia/packages/Plots/8K4be/src/backends.jl:31 [inlined]
  [7] plotlyjs()
    @ Plots ~/.julia/packages/Plots/8K4be/src/backends.jl:31
  [8] top-level scope
    @ In[4]:1
  [9] eval
    @ ./boot.jl:373 [inlined]
 [10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196

We can make a nice plot of the function by first choosing a much finer grid. However, the contour and surface plotting functions expect the transpose of mtx(\(f\)).

To emphasize departures from a zero level, use a colormap such as redsblues, and use clims to set balanced color differences.

m = 70;   x = range(0,2,length=m+1);
n = 50;   y = range(1,3,length=n+1);
F = [ f(x,y) for x in x, y in y ];

x y
x y f(x,y)

Parameterized surfaces

We are not limited to rectangles by tensor products. Many regions and surfaces may be parameterized by means of \(x(u,v)\), \(y(u,v)\), and \(z(u,v)\), where \(u\) and \(v\) lie in a rectangle. Such “logically rectangular” surfaces include the unit disk,

(13.1.5)\[\begin{split}\left\{ \begin{aligned} x &= u \cos v, \\ y &= u \sin v,\\ \end{aligned} \right. \qquad \qquad \left. \begin{aligned} 0 & \le u < 1, \\ 0 &\le v \le 2\pi, \end{aligned} \right.\end{split}\]

and the unit sphere,

(13.1.6)\[\begin{split}\left\{ \begin{aligned} x &= \cos u \sin v,\\ y &= \sin u \sin v,\\ z &= \cos v, \end{aligned} \right. \qquad \qquad \left. \begin{aligned} 0 & \le u < 2\pi, \\ 0 &\le v \le \pi. \end{aligned} \right.\end{split}\]
Demo 13.1.4

For a function given in polar form, such as \(f(r,\theta)=1-r^4\), construction of a function over the unit disk is straightforward using a grid in \((r,\theta)\) space.

r = range(0,1,length=41)
θ = range(0,2π,length=81)
F = [ 1-r^4 for r in r, θ in θ ]

    xlabel="r",ylabel="θ",title="A polar function")
r θ A polar function

Of course, we are used to seeing such plots over the \((x,y)\) plane, not the \((r,\theta)\) plane. For this we create matrices for the coordinate functions \(x\) and \(y\).

X = [ r*cos(θ) for r in r, θ in θ ]
Y = [ r*sin(θ) for r in r, θ in θ ]

    xlabel="x",ylabel="y",title="Function on the unit disk")
x y