Tensor-product discretizations - Fundamentals of Numerical Computation

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:

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

(13.1.1)

The × 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.

13.1.1Functions 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:

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

(13.1.4)

Example 13.1.1

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

\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}.

(13.1.5)

Example 13.1.2 (Functions on 2D grids)

Julia

MATLAB

Python

Example 13.1.2

Here is the grid from Example 13.1.1.

m = 4
x = range(0, 2, m+1)
n = 2
y = range(1, 3, 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]

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$ ).

using Plots
m, n = 80, 60
x = range(0, 2, m+1);
y = range(1, 3, n+1);
F = [f(x, y) for x in x, y in y]
contour(x, y, F';
    levels=21,  aspect_ratio=1,
    color=:redsblues,  clims=(-1, 1),
    xlabel="x",  ylabel="y" )

Example 13.1.2

Here is the grid from Example 13.1.1.

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

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

[mtx, X, Y] = tensorgrid(x, y);
f = @(x, y) cos(pi * x .* y - y);
F = mtx(f)

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$ ).

m = 80;  x = linspace(0, 2, m+1);
n = 60;  y = linspace(1, 3, n+1);
[mtx, X, Y] = tensorgrid(x, y);
F = mtx(f);

pcolor(X', Y', F')
shading interp
colormap(redsblues),  colorbar
axis equal
xlabel("x"),  ylabel("y")

Example 13.1.2

Here is the grid from Example 13.1.1.

m = 4
x = linspace(0, 2, m+1)
n = 2
y = linspace(1, 3, n+1)

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

f = lambda x, y: cos(pi * x * y - y)
F = array( [ [f(xi, yj) for yj in y] for xi in x ] )
print(F)

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$ ).

m, n = 80, 70
x = linspace(0, 2, m+1)
y = linspace(1, 3, n+1)
mtx, X, Y, _, _, _ = FNC.tensorgrid(x, y)
F = mtx(f)

pcolormesh(X.T, Y.T, F.T, cmap="RdBu", vmin=-1, vmax=1, shading="gouraud")
axis("equal"),  colorbar()
xlabel("$x$"),  ylabel("$y$");

13.1.2Parameterized 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,

\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.

(13.1.6)

and the unit sphere,

\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.

(13.1.7)

Example 13.1.3 (Functions on parameterized surfaces)

Julia

MATLAB

Python

Example 13.1.3

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, 41)
θ = range(0, 2π, 81)
F = [1 - r^4 for r in r, θ in θ]
plot(r, θ, F';
    legend=:none, 
    color=:viridis,  fill=true,
    xlabel="r",  ylabel="θ", 
    title="A polar function")

Of course, we are used to seeing such plots over the $(x,y)$ plane, not the $(r,\theta)$ plane.

In such functions the values along the line $r=0$ must be identical, and the values on the line $\theta=0$ should be identical to those on $\theta=2\pi$ . Otherwise the interpretation of the domain as the unit disk is nonsensical. If the function is defined in terms of $x$ and $y$ , then those can be defined in terms of $r$ and θ using (13.1.6).

Example 13.1.3

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 = linspace(0, 1, 41);
theta = linspace(0, 2*pi, 121);
[mtx, R, Theta] = tensorgrid(r, theta);
F = mtx(@(r, theta) 1 - r.^4);
clf,  colormap(parula)
contourf(R', Theta', F', 20)
shading flat
xlabel("r"),  ylabel("\theta"), 
title("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(Theta);  Y = R .* sin(Theta);
contourf(X', Y', F', 20)
axis equal,  shading interp  
xlabel('x'),  ylabel('y')
title('Function over the unit disk')

Example 13.1.3

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 = linspace(0, 1, 41)
theta = linspace(0, 2*pi, 121)
mtx, R, Theta, _, _, _ = FNC.tensorgrid(r, theta)

F = mtx(lambda r, theta: 1 - r**4)    

contourf(R.T, Theta.T, F.T, levels=20)
colorbar()
xlabel("$r$"),  ylabel("$\\theta$");

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, Y = R * cos(Theta), R * sin(Theta)
contourf(X.T, Y.T, F.T, levels=20)
colorbar(),  axis("equal")
xlabel("$x$"),  ylabel("$y$");

13.1.3Partial derivatives¶

In order to solve boundary-value problems in one dimension by collocation, we replaced an unknown function $u(x)$ by a vector of its values at selected nodes and discretized the derivatives in the equation using differentiation matrices. We use the same ideas in the 2D case: we represent a function by its values on a grid, and multiplication by differentiation matrices to construct discrete analogs of the partial derivatives $\frac{\partial u}{\partial x}$ and $\frac{\partial u}{\partial y}$ .

Consider first $\frac{\partial u}{\partial x}$ . In the definition of this partial derivative, the independent variable $y$ is held constant. Note that $y$ is constant within each column of $\mathbf{U} = \mtx(u)$ . Thus, we may regard a single column $\mathbf{u}_j$ as a discretized function of $x$ and, as usual, left-multiply by a differentiation matrix $\mathbf{D}_x$ such as (10.3.6). We need to do this for each column of $\mathbf{U}$ by $\mathbf{D}_x$ , which is accomplished by $\mathbf{D}_x \mathbf{U}$ . Altogether,

\mtx\left( \frac{\partial u}{\partial x} \right) \approx \mathbf{D}_x \, \mtx(u).

(13.1.8)

This relation is not an equality, because the left-hand side is a discretization of the exact partial derivative, while the right-hand side is a finite-difference approximation. Yet it is a natural analog for partial differentiation when we are given not $u(x,y)$ but only the grid value matrix $\mathbf{U}$ .

Now we tackle $\frac{\partial u}{\partial y}$ . Here the inactive coordinate $x$ is held fixed within each row of $\mathbf{U}$ . However, if we transpose $\mathbf{U}$ , then the roles of rows and columns are swapped, and now $y$ varies independently down each column. This is analogous to the situation for the $x$ -derivative, so we left-multiply by a finite-difference matrix $\mathbf{D}_y$ , and then transpose the entire result to restore the roles of $x$ and $y$ in the grid. Fortunately, linear algebra allows us to express the sequence transpose–left-multiply–transpose more compactly:

\mtx\left( \frac{\partial u}{\partial y} \right) \approx \Bigl(\mathbf{D}_y \mathbf{U}^T\Bigr)^T = \mtx(u)\, \mathbf{D}_y^T.

(13.1.9)

Keep in mind that the differentiation matrix $\mathbf{D}_x$ is based on the discretization $x_0,\ldots,x_m$ , and as such it must be $(m+1)\times (m+1)$ . On the other hand, $\mathbf{D}_y$ is based on $y_0,\ldots,y_n$ and is $(n+1)\times (n+1)$ . This is exactly what is needed dimensionally to make the products in (13.1.8) and (13.1.9) consistent. More subtly, if the differentiation is based on equispaced grids in each variable, the value of $h$ in a formula such as (5.4.8) will be different for $\mathbf{D}_x$ and $\mathbf{D}_y$ .

Example 13.1.4 (Numerical partial derivatives)

Julia

MATLAB

Python

Example 13.1.4

We define a function and, for reference, its two exact partial derivatives.

u = (x, y) -> sin(π * x * y - y);
∂u_∂x = (x, y) -> π * y * cos(πx * y - y);
∂u_∂y = (x, y) -> (π * x - 1) * cos(π * x * y - y);

We use an equispaced grid and second-order finite differences as implemented by diffmat2.

m = 80;
x, Dx, _ = FNC.diffmat2(m, [0, 2]);
n = 60;
y, Dy, _ = FNC.diffmat2(n, [1, 3]);
mtx = (f, x, y) -> [f(x, y) for x in x, y in y]
U = mtx(u, x, y)
∂xU = Dx * U
∂yU = U * Dy';

Now we compare the exact $\frac{\partial u}{\partial y}$ with its finite-difference approximation.

M = maximum(abs, ∂yU)    # find the range of the result
plot(layout=(1, 2), 
    aspect_ratio=1,   clims=(-M, M), 
    xlabel="x", ylabel="y")
contour!(x, y, mtx(∂u_∂y, x, y)';
    levels=15,  subplot=1,
    color=:redsblues,
    title="∂u/∂y")
contour!(x, y, ∂yU';
    levels=15,  subplot=2,
    color=:redsblues, 
    title="approximation")

To the eye there is little difference to be seen, though the results have no more than a few correct digits at these discretization sizes:

exact = mtx(∂u_∂y, x, y)
# Relative difference in Frobenius norm:
norm(exact - ∂yU) / norm(exact)

Example 13.1.4

We define a function and, for reference, its two exact partial derivatives.

u = @(x, y) sin(pi * x .* y - y);
du_dx = @(x, y) pi * y .* cos(pi * x .* y - y);
du_dy = @(x, y) (pi * x - 1) .* cos(pi * x .* y - y);

We will use an equispaced grid and second-order finite differences as implemented by diffmat2. First, we have a look at a plots of the exact partial derivatives.

m = 80;  [x, Dx] = diffmat2(m, [0, 2]);
n = 60;  [y, Dy] = diffmat2(n, [1, 4]);
[mtx, X, Y] = tensorgrid(x, y);
U = mtx(u);
dU_dX = mtx(du_dx);
dU_dY = mtx(du_dy);

clf,  subplot(1, 2, 1)
pcolor(X', Y', dU_dX')
axis equal,  shading interp
title('∂u/∂x')
subplot(1, 2, 2)
pcolor(X', Y', dU_dY')
axis equal,  shading interp
title('∂u/∂y')

Now we compare the exact partial derivatives with their finite-difference approximations. Since these are signed errors, we use a colormap that is symmetric around zero.

err = dU_dX - Dx * U;
subplot(1, 2, 1)
pcolor(X', Y', err')
colorbar,  clim([-.4, .4])
axis equal,  shading interp
title('error in ∂u/∂x')

err = dU_dY - U * Dy';
subplot(1,2,2)
pcolor(X', Y', err')
colorbar,  clim([-.1, .1])
axis equal,  shading interp
colormap(redsblues)
title('error in ∂u/∂y')

Not surprisingly, the errors are largest where the derivatives themselves are largest.

Example 13.1.4

We define a function and, for reference, its two exact partial derivatives.

u = lambda x, y: sin(pi * x * y - y)
du_dx = lambda x, y: pi * y * cos(pi * x * y - y)
du_dy = lambda x, y: (pi * x - 1) * cos(pi * x * y - y)

We will use an equispaced grid and second-order finite differences as implemented by diffmat2. First, we have a look at a plots of the exact partial derivatives.

m, n = 80, 60
x, Dx, Dxx = FNC.diffmat2(m, [0, 2])
y, Dy, Dyy = FNC.diffmat2(n, [1, 4])
mtx, X, Y, _, _, _ = FNC.tensorgrid(x, y)

U = mtx(u)
dU_dX = mtx(du_dx)
dU_dY = mtx(du_dy)

subplot(1, 2, 1)
contourf(X.T, Y.T, dU_dX.T)
title("$u$"),  axis("equal")
subplot(1, 2, 2)
contourf(X.T, Y.T, dU_dY.T)
title("$\\partial u/\\partial y$"),  axis("equal");

Now we compare the exact partial derivatives with their finite-difference approximations. Since these are signed errors, we use a colormap that is symmetric around zero.

subplot(1, 2, 1)
pcolormesh(X, Y, Dx @ U  - dU_dX, shading="gouraud", cmap="RdBu", vmin=-0.4, vmax=0.4)
colorbar()
title("error in $\\partial u/\\partial x$"),  axis("equal")
subplot(1, 2, 2)
pcolormesh(X, Y, U @ Dy.T - dU_dY, shading="gouraud", cmap="RdBu", vmin=-0.1, vmax=0.1)
colorbar()
title("error in $\\partial u/\\partial y$"),  axis("equal");

Not surprisingly, the errors are largest where the derivatives themselves are largest.

13.1.4Exercises¶

⌨ In each part, make side-by-side surface and contour plots of the given function over the given domain.
(a) $f(x,y) = 2y + e^{x-y}$ , $\quad[0,2]\times[-1,1]$
(b) $f(x,y) = \tanh[5(x+xy-y^3)]$ , $\quad[-2,2]\times[-1,1]$
(c) $f(x,y) = \exp \bigl[-6(x^2+y^2-1)^2 \bigr]$ , $\quad[-2,2]\times[-2,2]$
⌨ For each function in Exercise 1, make side-by-side surface plots of $f_x$ and $f_y$ using Chebyshev spectral differentiation.
⌨ For each function in Exercise 1, make a contour plot of the mixed derivative $f_{xy}$ using Chebyshev spectral differentiation.
⌨ In each case, make a plot of the function given in polar or Cartesian coordinates over the unit disk.
(a) $f(r,\theta) = r^2 - 2r\cos \theta$
(b) $f(r,\theta) = e^{-10r^2}$
(c) $f(x,y) = xy - 2 \sin (x)$
⌨ Plot $f(x,y,z)=x y - x z - y z$ as a function on the unit sphere.
⌨ Plot $f(x,y,z)=x y - x z - y z$ as a function on the cylinder $r=1$ for $-1\le z \le 2$ .

Preface

Two-dimensional problems

Preface

Two-dimensional diffusion and advection