Lecture 37: Time-evolution PDEs with periodic boundary conditions¶

This lecture investigates solving time-evolution PDEs with periodic boundary conditions. We will accomplish this initially using a semi-discretization: we discretize in space, and evolve in time exactly, using the matrix exponential.

We will consider two examples.

1. The transport equation $$u_t = u_\theta, u(0,\theta)=u_0(\theta)$$ which has exact solution: $u(t,\theta) = u_0(t+\theta)$. We will use this to compare our numerical approximation.

2. Heat equation $$u_t = u_{\theta\theta}, u(0,\theta)=u_0(\theta).$$

Discrete derivatives for periodic functions¶

We first construct discrete derivatives that incorporate boundary conditions. The trapezoidal approximation will now be 2π periodic. For example, consider the approximate $\cos(\cos(\theta))$ at $\theta_0,\ldots,\theta_n$:

In [76]:
n=5
θ=linspace(0,2π,n+1)
plot(θ,cos(cos(θ)));


Because $\cos\cos(\theta))$ is 2π periodic, we can extend it to the whole real line, as well as the piecewise affine approximation:

In [78]:
n=5
θ=linspace(0,2π,n+1)
plot(θ,cos(cos(θ)));
plot(θ-4π,cos(cos(θ-4π)));
plot(θ-2π,cos(cos(θ-2π)));
plot(θ+2π,cos(cos(θ+2π)));
plot(θ+4π,cos(cos(θ+2π)));


That way, we can also extend the derivative:

In [89]:
n=5
θ=linspace(0,2π,n+1)
h=1/n
dv=diff(cos(cos(θ)))/h
for k=1:length(dv)
plot([θ[k],θ[k+1]],[dv[k],dv[k]])
end

for k=length(dv)
plot([θ[k]-2π,θ[k+1]-2π],[dv[k],dv[k]])
end

for k=1
plot([θ[k]+2π,θ[k+1]+2π],[dv[k],dv[k]])
end


Thus if we use the left-hand discretization, that is, evaluate at $\theta_0+0,\ldots,\theta_{n-1}+0$, we wrap around. This gives us the $n \times n$ discrete derivative

$$D_n^L : {\hbox{Values at } \atop \theta_0,\ldots,\theta_{n-1}} \rightarrow {\hbox{Values at } \atop \theta_0,\ldots,\theta_{n-1}}$$

Defined by

$$D_n^L \triangleq {1 \over h} \begin{pmatrix} -1 & 1 \cr &-1 & 1 \cr &&\ddots & \ddots \cr &&&-1 & 1 \cr 1 &&&&-1\end{pmatrix}$$

We define this matrix as follows:

In [95]:
using PyPlot

function DL(n)
h=2π/n
ret=zeros(n,n)
for k=1:n-1
ret[k,k]=-1
ret[k,k+1]=1
end
ret[n,1]=1
ret[n,n]=-1
ret/h
end

n=1000
θ=linspace(0.,2π*(1-1/n),n)
norm(DL(n)*cos(θ)+sin(θ),Inf)

Out[95]:
0.003141582318193059

An equally valid construction is to evaluate at $\theta_0-0,\ldots,\theta_{n-1}-0$, in which case we get

$$D_n^R : {\hbox{Values at } \atop \theta_0,\ldots,\theta_{n-1}} \rightarrow {\hbox{Values at } \atop \theta_0,\ldots,\theta_{n-1}}$$

defined by

$$D_n^R \triangleq {1 \over h} \begin{pmatrix} 1 &&&& -1 \cr -1 & 1 \cr &-1 & 1 \cr &&\ddots & \ddots \cr &&&-1 & 1 \end{pmatrix}$$
In [97]:
function DR(n)
h=2π/n
ret=zeros(n,n)
for k=2:n
ret[k,k]=1
ret[k,k-1]=-1
end
ret[1,n]=-1
ret[1,1]=1
ret/h
end

n=1000
θ=linspace(0.,2π*(1-1/n),n)
norm(DR(n)*cos(θ)+sin(θ),Inf)

Out[97]:
0.003141582318193381

Solving PDEs with left-hand rule¶

Consider the transport equartion

$$u_t = u_\theta$$

We discretize the $\theta$-derivative using $D_n^L$, to obtain a time-dependent vector

$$\mathbf w(t) \approx \begin{pmatrix} u(t,\theta_0)\cr u(t,\theta_1) \cr \vdots \cr u(t,\theta_{n-1})\end{pmatrix}$$

that solves the ODE

$$\mathbf w'(t) = D_n^L \mathbf w(t)$$

with initial condition

$$\mathbf w(0) = \mathbf w_0 \triangleq \begin{pmatrix} u_0(\theta_0)\cr u_0(\theta_1) \cr \vdots \cr u_0(\theta_{n-1})\end{pmatrix}$$

The exact solution is given by the matrix exponential $e^{D_n^L t} \mathbf w_0$.

We can verify that this approximation converges:

In [14]:
u0=θ->cos(θ-0.9)

n=1000

θ=linspace(0.,2π*(1-1/n),n)
w0=u0(θ)
t=1.
w1=expm(DL(n)*t)*w0

plot(θ,w1)
plot(θ,u0(t+θ))

Out[14]:
1-element Array{Any,1}:
PyObject <matplotlib.lines.Line2D object at 0x30a189b90>

Solving PDEs with right-hand rule¶

We can also discretize the equation with $D_n^R:$

$$\mathbf w'(t) = D_n^R \mathbf w(t)$$

with initial condition

$$\mathbf w(0) = \mathbf w_0 \triangleq \begin{pmatrix} u_0(\theta_0)\cr u_0(\theta_1) \cr \vdots \cr u_0(\theta_{n-1})\end{pmatrix}$$

The exact solution is given by the matrix exponential $e^{D_n^R t} \mathbf w_0$. Unfortunately, it blows up and does not approximate the true solution:

In [34]:
u0=θ->cos(θ-0.9)

n=200
θ=linspace(0.,2π*(1-1/n),n)
w0=u0(θ)
t=1.
w1=expm(DR(n)*t)*w0
plot(θ,w1)
plot(θ,u0(t+θ))

Out[34]:
1-element Array{Any,1}:
PyObject <matplotlib.lines.Line2D object at 0x321381310>

Stability of time evolution¶

This phenomenon emphasizes the subtly of numerically solving time-evolution PDEs: getting a discretization "wrong" can have drastic effects. To understand what makes a "good" discretization, recall that, if $A=V \Lambda V^{-1}$, then

$$e^{At} = V e^{\Lambda t} V^{-1} = V \begin{pmatrix} e^{\lambda_1 t} \cr & \ddots \cr && e^{\lambda_n t} \end{pmatrix} V^{-1}.$$

Recall further that if $\lambda = q + i r$, then

$$e^{\lambda} = e^q e^{i r}.$$

Thus if $\lambda$ has positive real part, we get exponential increase and instability, and if $\lambda$ has negative real part, we have exponential decrease and stability.

Indeed, the eigenvalues of $D_n^L$ are all in the left-hand plane and so $e^{D_n^L t}$ is stable:

In [101]:
n=100
λ=eigvals(DL(n))
scatter(real(λ),imag(λ));


While the eigenvalues of $D_n^R$ are in the right-hand plane and so $e^{D_n^Rt} is unstable: In [103]: λ=eigvals(DR(n)) scatter(real(λ),imag(λ));  Excercise How would you discretize the PDE$u_t + u_{\theta} = 0$? Heat equation¶ Let's now consider the heat equation: $$u_t = u_{\theta\theta}$$ We have several choices for discretizing the second derivative: 1.$(D_n^L)^2$2.$(D_n^R)^2$3.$D_n^RD_n^L$4.$D_n^LD_n^R$Inspecting the eigenvalues, we see that 3 (and 4) are stable: In [105]: n=10 D2=DL(n)*DR(n) λ=eigvals(D2) scatter(real(λ),imag(λ));  Exercise Are 1 and 2 stable or unstable? We can therefore approximate the solution by $$e^{\tilde D_n^2 t} \mathbf w_0$$ where$\tilde D_n^2 = D_n^LD_n^R\$.

In [107]:
u0=θ->sin(θ-0.9)+1

n=200

D2=DL(n)*DR(n)

θ=linspace(0.,2π*(1-1/n),n)
w0=u0(θ)
t=1.
w1=expm(D2*t)*w0
t=2.
w2=expm(D2*t)*w0

plot(θ,w0)
plot(θ,w1)
plot(θ,w2);


If we investigate the derivatives closely, we observe that all eigenvalues are negative:

In [109]:
eigvals(D2)

Out[109]:
200-element Array{Float64,1}:
-4052.85
-4051.85
-4051.85
-4048.85
-4048.85
-4043.85
-4043.85
-4036.87
-4036.87
-4027.9
-4027.9
-4016.95
-4016.95
⋮
-35.8935
-24.9486
-24.9486
-15.979
-15.979
-8.99334
-8.99334
-3.99868
-3.99868
-0.999918
-0.999918
-2.52476e-14

Thus, eventually, we have artificial decay in the solution. To see this reliably, we have to implement expm ourselves using eigenvalue decomposition:

In [113]:
t=3000000000000.

λ,Q=eig(D2)
w=Q*diagm(exp(λ*t))*Q'*w0
plot(θ,w0)
plot(θ,w)   # should tend to 1, but some decay is introduced because   -2.52476e-14 < 0.;