6

I am trying to solve a linear equation system $\textbf{A}\textbf{x}=\textbf{b}$, e.g. a Poisson equation discretized in strong form, using biCGstab method.

Since there are only natural Neumann boundary conditions, at least one Dirichlet condition must be given. I simply add condition $\sum_{j} x=0$ to row $i$ in $\textbf{A}$. The result converges, but the error in the gradient field $\nabla x$ at $i$ is very large.

What is the reason for this, and how should I imposition this condition?

Kozuki
  • 99
  • 2
  • 5
  • 2
    Look at Bochev and Lehoucq, On the Finite Element Solution of the Pure Neumann Problem, SIAM Review 2005, http://epubs.siam.org/doi/abs/10.1137/S0036144503426074. – Christian Clason Jun 01 '15 at 14:05
  • 1
    Depending on your preconditioner, you may get away with not specifying such a condition. Krylov solvers will, in some sense, pick a pressure for you. – Bill Barth Jun 01 '15 at 14:31

3 Answers3

7

Pure Neumann problem is unique up to a constant. My two favourite solution strategies:

  1. Modifying the equation $-\Delta u = f$ to $-\Delta u + \varepsilon u = f$ for some small $\varepsilon>0$. If you perform reduced integration this corresponds roughly to adding $\varepsilon$ to the diagonal of your system matrix.

  2. Imposing $\int_\Omega u \,\mathrm{d}x=0$ using Lagrange multiplier. Then your system would be $$\begin{bmatrix} A & 1 \\ 1^T & 0 \end{bmatrix} \begin{bmatrix} x \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}.$$

knl
  • 2,096
  • 11
  • 17
  • Just a note: BiCGstab won't care, only CG requires symmetry and positive definiteness. (On the other hand, CG will work for the unmodified system, as long as the (discrete) system is consistent -- at least in exact arithmetic, and up to a constant.) – Christian Clason Jan 28 '16 at 19:33
  • Thanks for the clarification. I made a faulty assumption that this property is shared by CG and BiCGstab although I don't know the latter method well. – knl Jan 28 '16 at 20:12
  • @knl I tried to follow your approach in my last question, however it's not apparently working – FEGirl Oct 03 '21 at 10:52
  • The vector $1^{\text{T}}$ is only valid if each DOF has the same weight. A more general solution would be to consider the correct quadrature weights. Moreover, if the domain is not Cartesian one would also have to consider the correct Jacobian weights. Same has to be considered for method 1 using a diagonal "jitter". – ConvexHull Jul 27 '22 at 13:16
4

Including even one Dirichelt condition changes the problem you are trying to solve and will not give you the correct solution! You must fulfil the Discrete Compatibility Criteria, see e.g. first and second pages of: http://eprints.ma.man.ac.uk/894/2/covered/MIMS_ep2007_156_Sample_Chapter.pdf

It basically states that the summation of each element of b must be equal to zero for all-Neumann problems, otherwise your solution will drift as there is nothing holding it at the boundaries.

You simply need to add a line stating that b = b - mean(b) in your code before solving, if you are solving in double precision. In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria.

I attempted to solve an all-Neumann problem using CG by imposing a Dirichlet condition at a single point, as one is usually told to do anecdotally by e.g. some Prof., and then compared it to using a Discrete Cosine Transform approach. It does not yield the same result, using the Discrete Compatibility Criteria instead does.

DrHansGruber
  • 127
  • 5
  • Could you clarify what you mean by: "In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria."? – lightxbulb May 16 '21 at 15:06
0

I remember that the matrices from the Poisson equations with pure Neumann boundary conditions is singular. Is it right for your case ? If so, you cannot use the BiCGSTAB to solve the resulting linear systems.

Christian Clason
  • 12,301
  • 3
  • 48
  • 68
Hsien-Ming Ku
  • 129
  • 1
  • 9
  • Yes, it is singular, and that is why I have to add one Dirichlet bc. Because I only care about the gradient field. Maybe I can use QR solvers, but they are too slow. – Kozuki Jun 02 '15 at 04:59
  • Well, if the matrix is singular, maybe I suggest you use some regularization methods, such as Tikhonov regularization. – Hsien-Ming Ku Jun 02 '15 at 07:20
  • The resulting matrix will still be singular, even with the additional condition added. – shuhalo Jun 03 '15 at 07:44
  • @shuhalo but the new system converged fast, does that happen to a singular matrix? – Kozuki Jun 03 '15 at 13:48
  • @shuhalo If you specify one node as a dirichlet condition then I think the matrix will then be non-singular. e.g. if A =[1,-1;-1,1] and you make the first node dirichlet then now A = [1,0;-1,1] which is now non-singular. – James Jul 31 '15 at 23:06
  • 1
    I case someone still stumbles upon this answer, it is incorrect and Conjugate Gradients works for pure Neumann problems, see e.g. http://www.cs.sandia.gov/~pbboche/papers_pdf/2005SIREV.pdf – DrHansGruber Jan 21 '17 at 23:54
  • 1
    @DrHansGruber The link is dead. Could you tell me the title of the paper? Nvm, used the wayback machine, it's: ON THE FINITE ELEMENT SOLUTION OF THE PURE NEUMANN PROBLEM by Pavel Bochev. – lightxbulb May 15 '21 at 23:21