Why is the test function space in FEM chosen with homogeneous boundary conditions?

Question

It is so confusing, especially when I learns discontinuous galerkin method in broken Sobolev space and weak Dirichlet boundary condition.

If the trial function is chosen with homogeneous boundary conditions, how can I get a linear discrete equation corresponding to test basis functions with DOFs on the boundary? I mean if $v\in V_0 \subset V$, then the space of $u\in V$ would be "larger" than $V_0$. The final matrix would be rectangular rather than square. Extra conditions must be imposed to make it square and solvable.

My understanding now is that the homogeneous boundary condition is helpful in derivation of weak form from original differential equation, such as Poisson equation.

edit: When it comes to matrix, a $V_{0,h}\subset V_0$ and $V_{h} \subset V$ must be chosen first.

edit2: I am confused, sometimes the book may say bases of $V_0 $ are identical to $V,$ but the boundary condition is different.

Answer 1

This is more of a conceptual argument but may be helpful in understanding why the math works out this way. As a bonus, no matrix will be involved (because the problem is actually already with the weak formulation before any discretization).

What FEM does is to numerically approximate the solution $y\in V$ to some abstract operator equation $Ay=f$ for $A:V\to V^*$ for $f\in V^*$ and some abstract vector space $V$ and its dual $V^*$ (it's not really necessary for this argument to understand what a dual space is; it's enough to accept that it's the space for which the following works).

The first step is to see that $Ay-f = 0$ if and only if $\langle Ay-f,v\rangle = 0$ for all $v\in V$. (Intuitively, $0$ is the only vector that is orthogonal to every other vector.) Note that it's the same space $V$ for $y$ and $v$. It is crucial for this equivalence that $V$ is a vector space, in particular, that for $v_1,v_2 \in V$ you also have $v_1+v_2\in V$.

The next step is specific for PDEs. If $A$ is a differential operator such as the Laplacian $-\Delta$, it's easy to see that the solution $y$ is not unique unless you fix constants. The usual (but not the only) way is to fix the constant by fixing the value at the boundary. To include this in the above approach means only looking in the space of all functions that have the same values at the boundary (both for $y$ and for $v$, remember!) Now it should be fairly obvious that this doesn't work unless that value is $0$, because otherwise if you add two functions with the same boundary values, their boundary value will be different.

Of course, that's a too nasty restriction to take lying down. Luckily there's a trick: Pick a function $y_b$ that has the desired values at the boundary and satisfies $A y_b=0$, compute a solution $y_0$ to $Ay=f$ with zero boundary values as described above (including using test functions $v$ with zero boundary values), and set $y=y_0+y_b$. This is precisely what you are doing in computations, with some short cuts to avoid explicitly constructing $y_b$ first.

You mention the weak Dirichlet boundary conditions: Here you do not prescribe boundary conditions, but "penalize" non-zero values so that the solution will actually be (close to) zero. Since you do not have boundary conditions for the solution, you also shouldn't have boundary conditions for the test functions. (And for DG, there's additional difficulties since the weak form is an operator equation $A:V\to W^*$ for $W\neq V$, so the test functions have to come from a different space of functions that don't even need to have sensible values on the boundary.)

Answer 2

The important point, not mentioned in the answers above, is that the weak form is a variational formulation in which we ask the question: Can we find a function $u$ that minimizes the energy functional (or is a stationary point of some other functional) in such a way that any other potential solution $u+v$ does not minimize the energy (or is not a saddle point). The $v$ here becomes the test function in the weak formulation.

If we seek solutions among the set of functions that have boundary values $u|_\Gamma=g$, then clearly the perturbed functions also have to satisfy $(u+v)|_\Gamma=g$, which can only be the case if $v|_\Gamma=0$, i.e., satisfies homogeneous Dirichlet boundary conditions.

A lengthier discussion of this issue can be found in video lecture 21.5 and following here: http://www.math.colostate.edu/~bangerth/videos.html

Answer 3

While it is true that $V_0 \subset V$, these spaces are still both infinite dimensional, and there are no matrix sizes to speak of since no discretization has happened yet. To form matrices, you must choose finite dimensional subspaces $V_{0,h}$ and $V_h$ which "approximate" $V_0$ and $V$, in some sense. The size of the resulting matrices just depends on the dimension of $V_{0,h}$ and $V_h$. Therefore, if you use the same number of basis functions so that your finite-dimensional approximate spaces $V_{0,h}$ and $V_h$ have the same dimension, your matrices will end up square.

EDIT: Actually, the last sentence above is incorrect. You can get square matrices without having the exact same dimension for your subspaces; what matters is obviously the number of DOF you have, which also includes essential boundary conditions which are enforced in your problem. The following example should hopefully be more helpful than what I wrote above.

We seek a function $u(x)$ such that \begin{align} -u'' +u = f, \: x \in \Omega, \\ u'|_{\Gamma_N} = g_N, \: u|_{\Gamma_D} = g_D, \end{align} where $\Omega$ is some open interval whose boundary is partitioned $\partial \Omega = \Gamma_N \cup \Gamma_D$. One weak form for this problem is: find $u \in U$ such that \begin{align*} \int_{\Omega} v' u' + \int_{\Omega} vu = \int_{\Omega} vf - v|_{\Gamma_N} g_N , \: \forall v \in V, \end{align*} where \begin{align} V = \{ w \: | \: w \in H^1(\Omega), w|_{\Gamma_D} = 0 \}, \\ U = \{ w \: | \: w \in H^1(\Omega), w|_{\Gamma_D} = g_D \}. \end{align} Note that we set $w|_{\Gamma_D} = 0$ in the test function space $V$ to eliminate the $v|_{\Gamma_D} u''|_{\Gamma_D}$ term since no information about $u''$ is known on that part of the boundary. This does not affect the DOF boundary conditions which we will enforce later, since we can enforce them directly without using our test function space $V$.

Lets proceed by a Petrov-Galerkin method: we make finite-dimensional subspaces by specifying bases \begin{align*} U_m = \text{span} \{\phi_j\}_{j = 1}^m, \\ V_n = \text{span} \{\psi_i \}_{i = 1}^n. \end{align*} The weak problem reduces to problem of seeking some $u_m \in U_M$ for which the following $n+1$ equations hold: \begin{align*} \int_{\Omega} \psi_i' u_m' + \int_{\Omega} \psi_i u_m = \int_{\Omega} \psi_i f - \phi_i |_{\Gamma_N} g_N, \: i = 1, ... ,n \\ u_m|_{\Gamma_D} = g_D \text{ (enforce DOF bc directly without $V$) }. \end{align*} Expanding $u_m$ out in a basis, this becomes the linear system \begin{align*} \sum_{j = 1}^m \left( \int_{\Omega} \psi_i ' \phi_j ' + \psi_i \phi_j \right) u_j = \int_{\Omega} \psi_i f - \phi_i |_{\Gamma_N} g_N, \: i = 1, ..., n \\ \sum_{j = 1}^m u_j \phi_j = g_D. \end{align*} In matrix form (with m = 3, n = 2 for simplicity): \begin{align*} \begin{bmatrix} \int_{\Omega} (\psi_1' \phi_1' + \psi_1 \phi_1) & \int_{\Omega} (\psi_1' \phi_2' + \psi_1 \phi_2) & \int_{\Omega} (\psi_1' \phi_3' + \psi_1 \phi_3) \\ \int_{\Omega} (\psi_2' \phi_1' + \psi_2 \phi_1) & \int_{\Omega} (\psi_2' \phi_2' + \psi_2 \phi_2) & \int_{\Omega} (\psi_2' \phi_3' + \psi_2 \phi_3) \\ \phi_1 |_{\Gamma_D} & \phi_2 |_{\Gamma_D} & \phi_3 |_{\Gamma_D} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} \int_{\Omega} \psi_1 f - \phi_1|_{\Gamma_N} g_N \\ \int_{\Omega} \psi_2 f - \phi_2|_{\Gamma_N} g_N \\ g_D \end{bmatrix}. \end{align*} So, the last row to enforce essential boundary conditions extends the matrix to be square, as you mention. Once again, it was fine that we restricted our space $V$ to be zero on part of the boundary since the essential DOFs were enforced directly using the equation $\sum u_j \phi_j = g_D$, which does not involve any functions from $V$.