1

I have a question regarding the positive definiteness of the stiffness matrix. Specifically, I believe that it should be positive definite only when at least one Dirichlet point is given, so I would like a clarification on which is the precise moment that the stiffness matrix becomes positive definite. Let us consider a very simple model problem:

$$ -\Delta u(\pmb{x}) = f(\pmb{x}), x \in \Omega \\ u(\pmb{x}) = g(\pmb{x}), \pmb{x} \in \Gamma_D \\ \partial_n u(\pmb{x}) = 0, \pmb{x} \in \Gamma_N. $$

Using the divergence theorem this can be rewritten as:

$$ \int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega}fv \\ u(\pmb{x}) = g(\pmb{x}), \, \pmb{x} \in \Gamma_D. $$

A space $V_h \subset \mathcal{H}^1(\Omega)$ is given, with basis functions $\phi_1, \ldots, \phi_d, \phi_{d+1},\ldots,\phi_{d+m}$. Let $Dir = \{1,\ldots,d\}$ and $Ind = \{d+1,\ldots,d+m\}$. The solution would then be of the form:

$$u_h(\pmb{x}) = \sum_{k=1}^{d+m}u_k\phi_{k}(\pmb{x}), \pmb{x} \in \Omega,$$

and the linear system of equations that must be solved is:

$$ \sum_{k=1}^{d+m}u_k\int_{\Omega}\nabla\phi_i \cdot \nabla\phi_j = \int_{\Omega}f\phi_i, \, i \in Ind \\ u_i = g_i, \, i \in Dir. $$

The full stiffness matrix and right-hand side are:

$$ W_{ij} = \int_{\Omega}\nabla\phi_i \cdot \nabla\phi_j, \, i,j \in Dir\cup Ind \\ f_i = \int_{\Omega}f\phi, \, i \in Dir \cup Ind. $$

Then the system can be written more concisely as:

$$\pmb{W}|_{Ind \times Dir\cup Ind}\pmb{u} = \pmb{f}|_{Ind},$$

or if the Dirichlet columns are removed:

$$ \pmb{W}|_{Ind \times Ind}\pmb{u}|_{Ind} = \pmb{f}|_{Ind}-\pmb{W}|_{Ind \times Dir}\pmb{g}. $$

I know, that if even a single Dirichlet node is provided, then the system has a unique solution, which to me implies that $\pmb{W}|_{Ind \times Ind}$ is positive definite. However, if no Dirichlet node is provided, one is left with pure Neumann boundary conditions, and an additional constraint must be provided, for example (https://fenicsproject.org/olddocs/dolfin/2016.1.0/python/demo/documented/neumann-poisson/python/documentation.html):

$$\int_{\Omega}u = 0 \implies \sum_{k=1}^{m}u_k\int_{\Omega}\phi_k = 0$$.

To me, requiring such an additional constraint implies that the matrix $\pmb{W}|_{Ind \times Ind}$ is not positive definite and that it has a zero eigenvalue. Yet, the answer here: In FEM, why is the stiffness matrix positive definite? claims that it is positive definite without referring to the number of Dirichlet nodes. Clearly I am missing something, and I would like to understand what.

Edit:

At first glance $\pmb{W}|_{Ind \times Ind}$ seems oblivious to Dirichlet nodes, since it includes only terms involving basis functions that correspond to non-Dirichlet nodes. However, introducing Dirichlet nodes modifies the non-Dirichlet basis functions (even non-Dirichlet basis functions have to vanish at a Dirichlet node, thus implicitly affecting the terms in $\pmb{W}|_{Ind \times Ind}$). Is anyone aware of a reference that makes this argument rigorous? More precisely a proof that the matrix $\pmb{W}|_{Ind \times Ind}$ becomes positive-definite upon the introduction of such a Dirichlet point?

lightxbulb
  • 2,082
  • 8
  • 14
  • 2
    Short answer: You are right, in the absence of Dirichlet b.c.s Poisson problem is only positive semi-definite, and you have to add another condition to guarantee uniqueness of the solution. However, notice that, even if you don't add any other conditions, the solution is unique up to an additive constant and the nullspace only contains the unit element. – Abdullah Ali Sivas Feb 16 '21 at 23:37
  • @AbdullahAliSivas I am talking specifically about the matrix $\pmb{W}|_{Ind \times Ind}$. Is there a detailed formal proof of how introducing a Dirichlet node modifies this matrix in order to make it positive definite? I am missing that specific part. I also could not find at which point this detail is touched upon in: https://scicomp.stackexchange.com/questions/21423/in-fem-why-is-the-stiffness-matrix-positive-definite – lightxbulb Feb 17 '21 at 00:10
  • 1
    I guess I should have written the sentence "the nullspace only contains the unit element" differently. Using the notation in the accepted answer you linked, if there are no Dirichlet conditions, it is easy to show that $a(c,c)=0$ where $c$ is a non-zero constant function. Therefore, $a(u,v)$ is not coercive. You can similarly show that $\mathbf{1}^TW\mathbf{1}=0$, hence, $W$ is not positive definite in the absence of Dirichlet boundary conditions. – Abdullah Ali Sivas Feb 17 '21 at 01:15
  • Also in that answer, they choose their spaces carefully, namely $u,v\in H_{0}^1(\Omega)$ so Dirichlet boundary conditions are enforced implicitly -which is the standard procedure for most FE methods. – Abdullah Ali Sivas Feb 17 '21 at 01:17
  • @AbdullahAliSivas Is the point of your example to show a function $c$ which is not in $\mathcal{H}^1_{\Gamma_D}(\Omega)$ since it doesn't vanish on a Dirichlet boundary, and is instead only in $\mathcal{H}^1(\Omega)$. Then due to $V \subset \mathcal{H}^1$ it follows $\pmb{1}^T\pmb{W}\pmb{1} = \pmb{0}$ (if $c = 1$). I am still struggling to understand how this affects $\pmb{W}|{Ind \times Ind}$. Is it because by introducing Dirichlet nodes we modify the functions $\phi_i$ by requiring those to vanish at node not in $Ind$, so we get $\pmb{W}|{Ind \times Ind}$ to be positive definite from that? – lightxbulb Feb 17 '21 at 12:21
  • Not really. Consider the Poisson problem with only Neumann boundary conditions. Let $u$ solve that equation and $c\neq 0$ be a constant function. Then $u+c$ also solves that equation. Consider $a(c,c) = a(u+c-u,u+c-u) = a(u+c,u+c) - a(u,u) = 0$, which immediately shows us that $a(\cdot,\cdot)$ is not coercive, hence, the corresponding coefficient matrix is not positive definite (though it is positive semidefinite). By adding the extra condition $\int u =0$, you guarantee that $u+c$ is not solution, hence, $a(u,v)$ becomes coercive. Does that make sense? – Abdullah Ali Sivas Feb 17 '21 at 15:24
  • @AbdullahAliSivas I understand that the solution is defined up to a constant if no extra constraint is imposed. I am fine with that. I was more interested in the details. From what I have read, the bilinear form $a$ due to the Poisson equation applied to functions from $\mathcal{H}^1_{\Gamma_D}(\Omega)$ is in fact coercive. As I understand it, the functions $c$ simply isn't from $\mathcal{H}^1_{\Gamma_D}(\Omega)$ since it doesn't vanish on the Dirichlet boundary. That's my understanding as to why $a$ here is not coercive. Is that correct? – lightxbulb Feb 17 '21 at 15:34
  • 1
    Yes. $a(u,v)$ is not coercive on $H^1$, but it is coercive on $H^1_0$ (or in general, $H^1_g$ if $u=g$ on $\Gamma_D$) . – Abdullah Ali Sivas Feb 17 '21 at 21:30
  • @AbdullahAliSivas Thank you. Feel free to formulate this as an answer if you want so that I can accept it. – lightxbulb Feb 17 '21 at 21:36

0 Answers0