# Can a mixed boundary conditions in 1D Linear FE lead to non positive definite stiffness?

Consider the finite element discretisation of a general second order linear differential equation, with mixed boundary conditions at the end. The boundary conditions are given as $$\alpha_0u'(0)+\beta_0u(0)=\gamma_0$$ and $$\alpha_1u'(1)+\beta_1u(1)=\gamma_1$$. Where $$u$$ is the dependent variable and the domain is $$(0,1)$$. The reduced stiffness matrix after applying the boundary conditions also depends on the coefficients of the mixed boundary conditions: $$K_{11} \rightarrow K_{11}(\alpha_0,\beta_0)$$ and $$K_{nn} \rightarrow K_{nn}(\alpha_1,\beta_1)$$. Since $$\alpha_i,\beta_i$$ are arbitrary, can it so happen that for some values of these, the stiffness matrix is NOT positive definite?

• as far as I remember the Robin (and Neumann) BC's increase the diagonal value, discretize your BC and make use of it's value $K_{n,n+1}=f(K_{nn})$, where $f(K_{nn})$ is some function of the value inside your domain. So the matrix should remain pos def if it initially was. Commented Aug 28, 2023 at 15:04

Let's combine your $$\alpha,\beta,\gamma$$ into fewer symbols by diving your boundary conditions by $$\alpha$$.

For the general-dimensional case, $$-\Delta u = f, \qquad \text{in \Omega}, \\ \frac{\partial u}{\partial n} = \alpha u + g, \qquad \text{on \partial\Omega}, \\$$ the weak form you get after multiplying by a test function $$v$$, integrating over the domain integrating by parts, and using the boundary condition is $$a(u,v) = \int_\Omega \nabla u \cdot \nabla v - \int_{\partial \Omega}\alpha u v.$$ This form is for sure positive definite if $$\alpha<0$$. It is also positive definite if $$\alpha$$ is less than the Poincare constant on the domain, but if you make $$\alpha$$ large enough, then the whole thing becomes indefinite, and consequently your system matrix may also not be positive definite any more (even though it is symmetric).

Fortunately, physically correct models will always lead to cases where $$\alpha<0$$. For example, for the heat equation, $$\frac{\partial u}{\partial n}$$ is the outward heat flux from the object to a reservoir of temperature zero (assuming $$g=0$$), and for positive interior temperatures that means that there must be a positive outward heat flux -- that is, $$\alpha<0$$.

• Thank you. If I have a random differential equation which has a non PD stiffness, will the FEM still converge to the correct result? Commented Aug 28, 2023 at 23:17
• @aaquib This depends on your solver. For symmetric consistent positive semi-definite systems you can use conjugate gradients, for consistent symmetric indefinite systems you can use SYMMLQ, for symmetric indefinite inconsistent systems you can use MINRES, for non-symmetric consistent systems you can use BiCGStab, CGNE, for non-symmetric inconsistent systems you can use CGNR, LSQR. Inconsistent systems suggest that you wrote a PDE that has no solution (e.g. you specified contradictory boundary conditions), so that may be a mistake, however CGNR will find the projection with minimum norm. Commented Aug 29, 2023 at 8:00
• @aaquib There are of course many differential equations that, regardless of the boundary conditions, lead to indefinite linear systems. For many of these, the finite element method still works, but as lightxbulb states, you will need to use a different linear solver because Conjugate Gradients will no longer work for such systems. Commented Aug 29, 2023 at 16:04

You are right, positive definiteness is NOT always true. I checked it again for a simple case.

As a simple example, consider a 1-dimensional Laplace equation for the region $$0\leq x \leq 1$$: $$\begin{equation*} \frac{\partial^2 u(x)}{\partial x^2}=0. \end{equation*}$$ The Robin-like boundary conditions at $$x=0$$ and $$x=1$$ are specified in the question. The starting weighted form is as follows, $$\begin{equation*} \int_0^1\delta u(x)\frac{\partial^2 u(x)}{\partial x^2}dx=0, \end{equation*}$$ where $$\delta u(x)$$ is the arbitrary weighted function. As always, after the application of partial integration, we have the weak form: $$\begin{equation*} -\int_0^1\frac{\partial \delta u(x)}{\partial x}\frac{\partial u(x)}{\partial x}dx +\left.\delta u(x)\frac{\partial u(x)}{\partial x}\right|_{x=1} -\left.\delta u(x)\frac{\partial u(x)}{\partial x}\right|_{x=0} =0, \end{equation*}$$ whch is rewritten as, $$\begin{equation*} -\int_0^1\frac{\partial \delta u(x)}{\partial x}\frac{\partial u(x)}{\partial x}dx +\delta u(1)u'(1)-\delta u(0)u'(0)=0. \end{equation*}$$ Substitution of Robin-like B.C given in the question yields: $$\begin{equation*} \int_0^1\frac{\partial \delta u(x)}{\partial x}\frac{\partial u(x)}{\partial x}dx +\delta u(0)\frac{1}{\alpha_0}(-\beta_0u(0)+\gamma_0)-\delta u(1)\frac{1}{\alpha_1}(-\beta_1u(1)+\gamma_1)=0. \end{equation*}$$ If one divide $$[0,1]$$ region into $$(n-1)$$-finite elements, we get, \begin{align*} K_{11}\rightarrow K_{11}(\alpha_0,\beta_0)&=K_{11}+\frac{1}{\alpha_0}(-\beta_0), \\ K_{nn}\rightarrow K_{nn}(\alpha_1,\beta_1)&=K_{nn}-\frac{1}{\alpha_1}(-\beta_1), \\ K_{\text{others}}&\rightarrow\text{unchanged}. \end{align*} Thus the positive definiteness of K-matrix is NOT true under some boundary conditions.

• You should probably elaborate on how you reach the conclusion that $K$ is not positive definite from the above (e.g. trace being equal to the sum of the eigenvalues, and the eigenvalues of a positive definite matrix being positive). Commented Aug 28, 2023 at 12:42
• @lightxbulb Is it because you can always choose $\beta_0$ large enough? Commented Aug 28, 2023 at 13:09
• @aaquib You can make $-\beta_0/\alpha_0$ or/and $\beta_1/\alpha_1$ negative and large enough by absolute value to get $K_{11}$ or/and $K_{nn}$ negative. For positive definite matrices the diagonal elements have to be positive, so $K_{ii}<0$ would contradict this. Commented Aug 28, 2023 at 13:21
• Do each of the diagonal elements have to be positive? I thought it was the trace that mattered. Anyhow, does this affect in any way to convergence/stability of Robin BVPs? Because the only other non-positive definite problem I know is of buckling instability Commented Aug 28, 2023 at 13:34
• @aaquib if any diagonal $K_{ii}$ is negative, then choose the vector $v$ with $v_i=1$ and all other entries 0, then $v^\top K v = K_{ii} < 0$ so $K$ is not PD. Commented Aug 28, 2023 at 14:51