As far as I understand, you need to define boundary conditions in time and space to select a unique solution to a PDE and make it solvable. However, in ODEs I only need to specific the initial value in time which I think of as a boundary condition. If I now use the method of lines to transform a PDE into a system of ODEs, do I need the spatial boundary conditions? If yes, are they imposed by fixing some of the "lines" to given functions? If no, how can it be that they are first required but then somehow the need for them disappears?

I got this question because I did this with the 2D acoustic wave equation and the result looks reasonable, yet I did not impose any spatial boundary conditions. I have also read somewhere (that I sadly cannot find anymore) that there are implicit boundary conditions. Does that mean that me not specifying any conditions is equivalent to choosing some implicitly?


1 Answer 1


Let's take the example of the unsteady one-dimensional heat equation inside a solid on a domain $x\in[0,1]$: $$\partial_t u - D\partial_{xx} u = 0$$ with the initial profil $u(0,x) = u_0(x)$ at $t=0$, and Dirchlet boundary conditions enforcing that the wall at $x=0$ (respectively $x=1$) is at temperature $u_{L}$ (respectively $u_{R}$). If we use discretise our solid with a uniformely distributed set of points $i \in [1,N]$ at positions $x_i = i \Delta x$ with $\Delta x = 1/(N+1)$ and use second order centered finite differences for the diffusion term, we obtain the following semi-discrete set of ODEs: $$d_t u_i = D\frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}$$ for $i \in [2, N-1]$. The boundary points $i=1$ and $i=N$ must be taken care of separately to ensure the boundary conditions are respected. Let's take the example of the boundary at $x=0$. The temperature $u(x=0)$ is enforced to be 0. Our first point is as $x=1/(n+1)$, therefore we can write: $$d_t u_1 = D\frac{u_{2} - 2u_1 + u_{L}}{\Delta x^2}$$

Similarly for our last point: $$d_t u_N = D\frac{u_{N-1} - 2u_N + u_{R}}{\Delta x^2}$$

We have to use these formulations, otherwise we would not even be able to compute the centered finite difference approximation of the diffusion terms at these outer points. An error would for example be to aproximate these terms with a decentered fintie difference formula, e.g for the point $i=1$: $$ d_t u_1 = D\frac{ 2u_1 - 5u_2 + 4 u_3 - u_4}{\Delta x^2} $$ which uses a second-order forward finite difference expression of the diffusion term. But then, there is not guarantee that the temperature at $x=0$ will remain at $u_L$, i.e. your system will most likely have an unphysical behaviour.

Sometimes we rather use a Neumann boudary condition to specify a given heat flux at the boundary. Let's take the case of a Neumann BC at $x=0$, with a zero heat flux, i.e. $\frac{\partial u}{\partial x}(x=0) = 0$. This condition can be discretised as $\frac{u_1-u_L}{\Delta x}=0$, which means that $u_L(t) = u_1(t)$ at all times $t$. You can then inject this value of $u_L$ in the second equation of this answer.

Technically, a discretized PDE is integrated in time as follows:

  • the time integration routines gives you the discrete solution profile $u_i, i \in [1,N]$ at time $t_n$
  • we compute the boundary values $u_L$ and $u_R$ so as to enforce the proper boundary conditions
  • we compute the time derivatives of our discrete solution values and return them to time integration routine.

I don't know how you've managed to integrate your system without specifying BCs at all times, maybe its just luck in your formulation or the particular equation you consider. But in the case of the discrete heat equation I've treated here, I see 3 potential ways of having this problem:

  • you forget to specify $d_t u_1$ and then you actually enforce a Dirichlet boundary condition that reads $u(t,x_1)= u_0(x_1)$, i.e. the wall temperature remains equal to its initial value

  • you forget the term $u_L$ in the third equation, which the means you implicitly assume $u_L=0$ .

  • you used the forward finite difference expression, thus omitting the BC value. In this case, we can see what the implcitly asusmed of $u_L$ is by equating the 3rd and 5th equation and we obtain: $u_L = -4u_1 + 6u_2 - 4u_3 + u_4$, that is $u_L$ is relatively "unpredictable", and consequently so will be $u_1$, and $u_2$ and so on...

  • $\begingroup$ Thank you very much for the detailed answer! One thing I did not mention in the question is that I discretized with finite elements in space with triangular elements. Then I get two the mass and stiffness matrices and the equation $Au_{tt} = Bu^{(t)}$ that can be solved for the second t-derivates which I then give to the ODE solver. And here I just used the $A$ and $B$ that came out of the computation. However, my suspicion is that the entries for the boundary points are wrong because they are not surrounded by triangles on all sides. They are on the boundary after all. $\endgroup$
    – Marten
    Nov 23, 2020 at 12:56
  • $\begingroup$ I haven't done any finite elements for years, so I'd rather not answer on that :/ As far as I remember the boundary conditions are already taken into account when you write your weak formulation of the problem. $\endgroup$
    – Laurent90
    Nov 23, 2020 at 12:58
  • 3
    $\begingroup$ From reading this sfepy.org/doc-devel/solving_pdes_by_fem.html, it looks like if you don't do anything for implementing BC in FEM you end up with Neumann zero flux BC (what they call "natural" BC); at least that's what it looks like in the example discussed there, not sure how general this is. $\endgroup$ Nov 23, 2020 at 19:24

Your Answer

