Two-dimensional diffusion and advection

We next describe how to apply the method of lines to PDEs of the form

The PDE may be of either parabolic or hyperbolic type, with the primary difference being potential restrictions on the time step size. To keep descriptions and implementations relatively simple, we will consider only periodic conditions, or Dirichlet boundary conditions imposed on all the edges of the rectangular domain.

As described in Tensor-product discretizations, the rectangular domain is discretized by a grid $(x_i,y_j)$ for $i=0,\ldots,m$ and $j=0,\ldots,n$. The solution is semidiscretized as a matrix $\mathbf{U}(t)$ such that $U_{ij}\approx u(x_i,y_j,t)$. Terms involving the spatial derivatives of $u$ are readily replaced by discrete counterparts: $\mathbf{D}_x\mathbf{U}$ for $u_x$, $\mathbf{U}\mathbf{D}_{y}^T$ for $u_y$, and so on.

13.2.1Matrix and vector shapes¶

Our destination is an IVP that can be solved by a Runge–Kutta or multistep solver. These solvers are intended for vector problems, but our unknowns naturally have a matrix shape, which is the most convenient for the differentiation formulas (13.1.8) and (13.1.9). Fortunately, it’s easy to translate back and forth between a matrix and an equivalent vector.

The function vec is built-in to Julia, whereas unvec is a particular use case of the reshape function.

13.2.2Periodic end conditions¶

If the boundary conditions are periodic, then the unknowns in the method of lines are the elements of the matrix $\mathbf{U}(t)$ representing grid values of the numerical solution. For the purposes of an IVP solution, this matrix is equivalent to the vector $\mathbf{u}(t)$ defined as $\mathbf{u}=\operatorname{vec}(\mathbf{U})$.

13.2.3Dirichlet conditions¶

In Boundaries we coped with boundary conditions by removing the boundary values from the vector of unknowns being solved in the semidiscretized ODE. Each evaluation of the time derivative required us to extend the values to include the boundaries before applying differentiation matrices in space. We proceed similarly here, except that we have changes in the shape as well as boundary conditions to consider.

Suppose we are given a matrix $\mathbf{U}$ that represents the solution on an $(m+1)\times (n+1)$ grid, including boundary values. Then we define

where

is $(m-1)\times (m+1)$, and $\mathbf{E}_y$ is analogous but of size $(n-1)\times (n+1)$. The left multiplication in (13.2.4) deletes the first and last row of $\mathbf{U}$ and the right multiplication deletes its first and last column. All that remains, then, are the interior values, which are converted into a vector by the vec operator.

For the inverse transformation, let us be given a vector $\mathbf{w}$ of interior solution values. Then we define

This operator reshapes the vector to a grid of interior values, then appends one extra zero row and column on each side of the grid.^[1]

Now suppose the ODE unknowns for the interior solution values are in the vector $\mathbf{w}(t)$. When we form $\operatorname{unpack}(\mathbf{w})$, we reinterpret the values on the tensor-product grid and then extend these values to zero around the boundary. If the boundary values are given as $g(x,y)$, then $g$ has to be evaluated at the boundary nodes of the grid and inserted into the grid function matrix. Then the grid values are used to compute partial derivatives in $x$ and $y$, the discrete form of the PDE is evaluated, and we pack the result as the computation of $\mathbf{w}'$.

The wave equation introduces a little additional complexity. First, we write the 2D wave equation $u_{tt}=c^2(u_{xx}+u_{yy})$ in first-order form as

Now the grid unknowns are a pair of matrices $\mathbf{U}(t)$ and $\mathbf{V}(t)$. Typical boundary conditions would prescribe $u$ on all of the boundary and let $v$ be unspecified. Since the boundary values of $\mathbf{U}$ are prescribed, those values are omitted from the semidiscretization IVP, while all of $\mathbf{V}$ is included. All of these unknowns need to be packed into and unpacked from a single vector $\mathbf{w}(t)$ for the IVP solver.

13.2.4Exercises¶

1. ⌨ For the given $u(x,y)$, make a plot of the given quantity on the square $[-2,2]^2$ using appropriate differentiation matrices.

   (a) $u(x,y) = \exp(x-y^2)$; plot $u_{xx}+u_{yy}$

   (b) $u(x,y) =\cos (\pi x)+\sin (\pi y)$; plot $u_x+u_y$

   (c) $u(x,y) =\exp(-x^2-4y^2)$; plot $x u_y$

2. ⌨ Following Demo 13.2.2 as a model, solve the Allen–Cahn equation $u_t=u(1-u^2)+0.001(u_{xx}+u_{yy})$ on the square $[-1,1]^2$ with periodic conditions, taking $u(x,y,0)=\sin(\pi x)\cos(2\pi y)$. Use $m=n=60$ to solve up to $t=4$, and make an animation of the result.

3. ⌨ Following Demo 13.2.3 as a model, solve $u_t=y u_x-u_y+0.03(u_{xx}+u_{yy})$ on the square $[-1,1]^2$, with $u(x,y,0)=(1-x^2)(1-y^2)$ and homogeneous Dirichlet boundary conditions. Use $m=n=40$ to solve up to $t=2$, and make an animation of the result.

4. ⌨ Following Demo 13.2.4 as a model, solve $u_{tt}=u_{xx}+u_{yy}+\cos(7t)$ on the square $[-1,1]^2$, with $u(x,y,0)=x(1-x^6)(1-y^2)$, $u_t(x,y,0)=0$, subject to homogeneous Dirichlet boundary conditions. Take $m=n=60$ to solve between $t=0$ and $t=12$, and make an animation of the result.

5. From Maxwell’s equations we can find a way to convert the wave equation to a first-order form that, unlike (13.2.7), uses only first-order derivatives in space:

   $$\begin{split} u_t &= c^2(v_y - w_x),\\ v_t &= u_y, \\ w_t &= -u_x, \end{split}$$

   (13.2.8)

   subject to $u=0$ on the boundary.

   (a) ✍ Show that a solution of (13.2.8) satisfies $u_t=c^2(u_{xx}+u_{yy})$.

   (b) ⌨ Solve (13.2.8) with $c=2$ in the rectangle $x\in[-3,3]$, $y\in[-1,1]$, $u(x,y,0) = \exp(x-x^2)(9-x^2)(1-y^2)$, and $v=w=0$ at $t=0$. Use $m=50$ for $x$ and $n=25$ for $y$, solve for $0\le t \le 6$, and make an animation of the solution.