```
```{raw} latex
%%start demo%%
```
```{code-cell}
m = 2; n = 3;
V = rand(1:9,m,n);
v = vec(V)
```
The `unvec` operation is the inverse of vec.
```{code-cell}
unvec = z -> reshape(z,m,n)
unvec(v)
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Periodic 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})$.
```{index} heat equation
```
(demo-diffadv-heat)=
```{proof:demo}
```
```{raw} html
```
```{raw} latex
%%start demo%%
```
We will solve a 2D heat equation, $u_t = 0.1(u_{xx} + u_{yy})$, on the square $[-1,1]\times[-1,1]$, with periodic behavior in both directions.
```{code-cell}
m = 60; x,Dx,Dxx = FNC.diffper(m,[-1,1]);
n = 25; y,Dy,Dyy = FNC.diffper(n,[-1,1]);
mtx = f -> [ f(x,y) for x in x, y in y ]
unvec = z -> reshape(z,m,n);
```
Note that the initial condition should also be periodic on the domain.
```{code-cell}
u_init = (x,y)->sin(4*π*x)*exp(cos(π*y))
U₀ = mtx(u_init)
M = maximum(abs,U₀)
contour(x,y,U₀',color=:redsblues,clims=(-M,M),aspect_ratio=1,
xaxis=("x",(-1,1)),yaxis=("y",(-1,1)),title="Initial condition")
```
This function computes the time derivative for the unknowns. The actual calculations take place using the matrix shape.
```{code-cell}
function dudt(u,α,t)
U = unvec(u)
Uxx = Dxx*U; Uyy = U*Dyy' # 2nd partials
dUdt = α*(Uxx + Uyy); # PDE
return vec(dUdt)
end;
```
Since this problem is parabolic, a stiff integrator is appropriate.
```{code-cell}
IVP = ODEProblem(dudt,vec(U₀),(0,0.2),0.1)
sol = solve(IVP,Rodas4P());
```
::::{panels}
:column: col-7 left-side
:card: border-0 shadow-none
```{raw} latex
\begin{minipage}[t]{0.5\textwidth}
```
Here is an animation of the solution.
```{raw} latex
\end{minipage}\hfill
```
---
:column: col-5 right-side
:card: shadow-none comment
```{raw} latex
\begin{minipage}[t]{0.4\textwidth}\begin{mdframed}[default]\small
```
Here `clims` are set so that colors remain at fixed values throughout the animation.
```{raw} latex
\end{mdframed}\end{minipage}
```
::::
```{code-cell}
:tags: [hide-input]
anim = @animate for t in range(0,0.2,length=81)
surface(x,y,unvec(sol(t))',color=:redsblues,clims=(-M,M),
xaxis=(L"x",(-1,1)),yaxis=(L"y",(-1,1)),zlims=(-M,M),
title=@sprintf("Heat equation, t=%.3f",t),
dpi=100,colorbar=:none)
end
closeall(); mp4(anim,"figures/diffadv-heat.mp4")
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Dirichlet conditions
```{index} boundary conditions; numerical implementation of
```
In {numref}`section-diffusion-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
:::{math}
:label: tpdeletion
\operatorname{pack}(\mathbf{U}) = \operatorname{vec}(\mathbf{E}_x \mathbf{U} \mathbf{E}_y^T),
:::
where
$$
\mathbf{E}_x = \begin{bmatrix}
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 \\
& & & \ddots & & \\
0 & 0 & 0 & \cdots & 1 & 0
\end{bmatrix}
$$
is $(m-1)\times (m+1)$, and $\mathbf{E}_y$ is analogous but of size $(n-1)\times (n+1)$. The left multiplication in {eq}`tpdeletion` 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
$$
\operatorname{unpack}(\mathbf{w}) = \mathbf{E}_x^T \cdot \operatorname{unvec}(\mathbf{w}) \cdot \mathbf{E}_y.
$$
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.[^jacobian]
[^jacobian]: You might wonder why we use linear algebra to define the extension and deletion of boundary values rather than directly accessing row and column indices in the grid function. The linear algebra approach allows `DifferentialEquations` to compute the Jacobian matrix of the implicit IVP solver quickly using *automatic differentiation* tools, greatly speeding up the solution process. Since the matrices in our expressions are sparse, multiplications by them do not affect running time perceptibly.
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}'$.
```{index} advection-diffusion equation
```
(demo-diffadv-advdiff)=
```{proof:demo}
```
```{raw} html
```
```{raw} latex
%%start demo%%
```
We will solve an advection-diffusion problem, $u_t + u_x = 1 + \epsilon(u_{xx} + u_{yy})$, where $u=0$ on the boundary of the square $[-1,1]^2$. The outline of our approach is based on {numref}`Function {number} ` for parabolic PDEs in one space dimension.
The first step is to define a discretization of the domain. We use the result to discretize an initial condition.
```{code-cell}
m = 50; n = 36;
x,Dx,Dxx = FNC.diffcheb(m,[-1,1])
y,Dy,Dyy = FNC.diffcheb(n,[-1,1])
U₀ = [ (1+y)*(1-x)^4*(1+x)^2*(1-y^4) for x in x, y in y ];
```
Next, we define the pack and unpack functions for the grid unknowns.
```{code-cell}
Ex = [zeros(m-1) I(m-1) zeros(m-1)]
Ey = [zeros(n-1) I(n-1) zeros(n-1)]
unvec = u -> reshape(u,m-1,n-1)
pack = U -> vec(Ex*U*Ey')
unpack = w -> Ex'*unvec(w)*Ey;
```
Now we can define and solve the IVP using a stiff solver.
```{code-cell}
function dwdt(w,ϵ,t)
U = unpack(w)
Ux,Uxx = Dx*U , Dxx*U
Uyy = U*Dyy'
dUdt = @. 1 - Ux + ϵ*(Uxx + Uyy)
return pack(dUdt)
end
IVP = ODEProblem(dwdt,pack(U₀),(0.,2),0.05)
w = solve(IVP,Rodas4P());
```
When we evaluate `w` at a particular value of $t$, we get a vector of the interior grid values. The same `unpack` function above converts this to a complete matrix of grid values.
```{code-cell}
anim = @animate for t in 0:0.02:2
U = unpack(w(t))
surface(x,y,U',layout=(1,2),size=(640,320),
xlabel=L"x",ylabel=L"y",zaxis=((0,2),L"u(x,y)"),
color=:viridis,alpha=0.66,clims=(0,2),colorbar=:none,
title="Advection-diffusion",dpi=100 )
contour!(x,y,U',levels=24,aspect_ratio=1,
subplot=2,xlabel=L"x",ylabel=L"y",
color=:viridis,clims=(0,2),colorbar=:none,
title=@sprintf("t = %.2f",t) )
end
closeall(); mp4(anim,"figures/diffadv-advdiff.mp4")
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
```{index} wave equation
```
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
:::{math}
:label: wave2dfirst
u_t &= v, \\
v_t &= c^2(u_{xx}+u_{yy}).
:::
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.
(demo-diffadv-wave)=
```{proof:demo}
```
```{raw} html
```
```{raw} latex
%%start demo%%
```
We solve the wave equation with $c=1$ on the square $[-2,2]\times[-2,2]$, where $u=0$ on the boundary. We start with the discretization and initial condition.
```{code-cell}
m = 40; n = 40;
x,Dx,Dxx = FNC.diffcheb(m,[-2,2])
y,Dy,Dyy = FNC.diffcheb(n,[-2,2])
U₀ = [ (x+0.2)*exp(-12*(x^2+y^2)) for x in x, y in y ]
V₀ = zeros(size(U₀));
```
We again define two functions for transforming between the matrix and vector representations. Unpacking extends the boundary values of $u$ by zeros, which would need to be modified if the PDE had a nonhomogeneous condition. Also note that there are two different sizes of unvec operations.
```{code-cell}
Ex = [zeros(m-1) I(m-1) zeros(m-1)]
Ey = [zeros(n-1) I(n-1) zeros(n-1)]
unvec_u = u -> reshape(u,m-1,n-1)
unvec_v = u -> reshape(u,m+1,n+1)
N = (m-1)*(n-1)
pack = (U,V) -> [ vec(Ex*U*Ey'); vec(V) ]
unpack = w -> ( Ex'*unvec_u(w[1:N])*Ey, unvec_v(w[N+1:end]) );
```
We can now define and solve the IVP. Since this problem is hyperbolic, not parabolic, a nonstiff integrator is faster than a stiff one.
```{code-cell}
function dwdt(w,c,t)
U,V = unpack(w)
dUdt = V
dVdt = c^2*(Dxx*U + U*Dyy')
return pack(dUdt,dVdt)
end
IVP = ODEProblem(dwdt,pack(U₀,V₀),(0,4.0),1)
sol = solve(IVP,Tsit5());
```
```{code-cell}
anim = @animate for t in 0:4/100:4
U,_ = unpack(sol(t))
surface(x,y,U',layout=(1,2),size=(640,320),
xlabel=L"x",ylabel=L"y",zaxis=((-0.1,0.1),L"u(x,y)"),
color=:redsblues,alpha=0.66,clims=(-0.1,0.1),colorbar=:none,
title="Wave equation",dpi=100 )
contour!(x,y,U',levels=24,aspect_ratio=1,
subplot=2,xlabel=L"x",ylabel=L"y",
color=:redsblues,clims=(-0.1,0.1),colorbar=:none,
title=@sprintf("t = %.2f",t) )
end
closeall(); mp4(anim,"figures/diffadv-wave.mp4")
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Exercises
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 {numref}`Demo %s