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

Question

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.

Answer 1

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.

Answer 2

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.