How do you apply boundary conditions in a time-stepping problem?

It looks to me like a very common problem, yet I haven't been able to find any practical guide on the subject despite many hours searching.

Here is a clearer statement of my question:

I have a system of (linear) PDEs in the variables $x$ and $t$. I can use the method of lines to discretise the space coordinate $x$ which yields a system of ODEs in time: $$\frac{d\vec{u}}{dt}=A\vec{u}+\vec{b}$$ $A$ is a square matrix.

Now, to me, the straightforward way to solve the above on a computer is to use an explicit time-stepping method of some sort. So I write:

$$\vec{u}_{n+1}=\vec{u}_{n}+\Delta t(A\vec{u}_{n}+\vec{b}_n)$$

Now, how do I enforce boundary conditions in all this ?

In their most general form, these are $$B\vec{u}=0$$ $B$ is a rectangular matrix.

I don't really see how this can fit into the above. Any help would be super appreciated.

• Boundary conditions are implicitly defined in the discrete spatial derivative operator $A$ – HBR Dec 20 '17 at 22:36
• How does it work in practice ? – jrekier Dec 20 '17 at 22:49
• This is discussed in every book on numerical methods for PDEs. – David Ketcheson Dec 21 '17 at 4:37
• David, your comment is really not helping. I have browsed through a fair number of such books and haven't found an explanation I could understand. Hence the question. For your comment to be of any value, you should at least provide a reference for one such book. – jrekier Dec 21 '17 at 8:39
• @jerkier a good introductory reference to numerical PDEs is this book by Randall Leveque. you could start by reading chapter 2.4 that deals with Dirichlet boundary conditions for a simple steady-state problem (in your notation, that'd be $Ax = b$), then work your way through the more advanced examples. – GoHokies Dec 21 '17 at 12:52

As HBR mentioned, the boundary conditions can often be immediately incorporated into $A$ and $b$. For example, suppose we wish to solve the 1D heat equation with Dirchilet boundary conditions

$$u_t = \sigma u_{xx} + h(x,t), \quad u(a,t) = f(t), \quad u(b,t) = g(t).$$

We then discretize $u_j^n \approx u(x_j, t_n)$ where $x_j = a+j \Delta x$ and $t_n = n\Delta t$. Here $j = 0,1,2,\ldots,N+1$ and $\Delta x = (b-a)/(N+1)$. Note that by our boundary conditions $u_0^n = f(t_n)$ and $u_{N+1}^n = g(t_n)$. As these values are completely determined for all time, we don't even include them in our linear system. For interior gridpoints ($2 \le j \le N-1$), we take a second-order centered difference

$$u_{xx}(x_j,t_n) \approx \frac{u_{j-1}^n - 2u_j^n + u_{j+1}^n}{\Delta x^2}. \tag{1}$$

For the boundary-adjacent gridpoints, we get

\begin{align} u_{xx}(x_1,t_n) &\approx \frac{u_{0}^n - 2u_1^n + u_{2}^n}{\Delta x^2}=\frac{f(t_n) - 2u_1^n + u_{2}^n}{\Delta x^2}, \tag{2} \\ u_{xx}(x_N,t_n) &\approx \frac{u_{N-1}^n - 2u_{N}^n + u_{N+1}^n}{\Delta x^2}=\frac{u_{N-1}^n - 2u_{N}^n + g(t_n)}{\Delta x^2}. \tag{3} \end{align}

Collecting our discretization, we get

$$\frac{d}{dt} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{n-1} \\ u_n \end{pmatrix} = \frac{1}{\Delta x^2}\begin{pmatrix} -2 & 1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & \cdots & 0 & 0& 0 \\ 0 & 1 & -2 & 1 & \cdots & 0& 0 & 0 \\ 0 & 0 & 1 & -2 & \cdots & 0 & 0& 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \cdots & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 & -2 \end{pmatrix}\begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ \vdots \\ u_{n-1} \\ u_n \end{pmatrix} + \begin{pmatrix} f(t)/\Delta x^2 + h(x_1,t) \\ h(x_2,t) \\ h(x_3,t) \\ h(x_4,t) \\ \vdots \\ h(x_{n-1},t) \\ g(t)/\Delta x^2 + h(x_n, t) \end{pmatrix}. \tag{4}$$

As you can see, we have baked our boundary conditions into the matrix $A$ and the vector $b$. (Make sure you can see how we got from the (1)-(3) to (4).) Similar approaches work for other PDEs and other discretization. Handling boundary conditions for hyperbolic problems is more tricky though, and you should always try to be very careful when treating the boundary in your discretization.

• Hi @eepperly, and thanks for your detailed answer. Much appreciated. It seems to me like this wouldn't be directly applicable with methods based on spectral differentiation (spectral methods). Am I right ? I will try your suggestion using collocation points though. Cheers ! – jrekier Dec 21 '17 at 8:49
• @jerkier a good reference for spectral methods is Jan Hesthaven's book. I can't point you right now to the specific section dealing with BCs (don't have the book at hand), but could do so later today if you're interested. – GoHokies Dec 21 '17 at 13:00

A simple way which should work well for an implicit time-integrator that can handle stiff equations is to write your BC's as relaxation ODEs. For example, if we want to enforce $u_j=u_{j0}$ where $u_j$ is the value of the unknown at a boundary grid point $j$, then use $\frac{d}{dt} u_j = - \lambda_{rlx} (u_j - u_{j0})$, where $\lambda_{rlx}$ is a "relaxation rate" chosen to be fast compared to time scales in the operator A.

The advantage of this method is it is very easy to implement. This way your BC is maintained not exactly (depending on the chosen $\lambda_{rlx}$), but the error can be made small, and sometimes it is actually desired to have some freedom in the boundary values, e.g., if your main goal is achieving a steady state.

• I will give it a try right away. Thx! – jrekier Dec 21 '17 at 8:50