1
$\begingroup$

The question is in the context of iterative numerical solution of large PDE systems with Finite Differences or Finite Elements:

Stating the Poisson equation with Neumann boundary conditions will lead to a singular system because it is invariant when adding a constant function. An often proposed option is to just "pin" the system to a fixed value in one point, which leads to a well-defined system.

In my case, I'm rather interested in a solution where the constant part vanishes, so "pinning" a certain point does not help.

Furthermore, I'd like to use an iterative solver, and those seem to be significantly unhappy with a singular system. In the PETSc documentation (I may or may not use PETSc then) I found the option to set a null space of the system.

So how should I deal with that singluar system? I'm interested in some theoretical views as well as practical solutions using linear solver libraries.

$\endgroup$
2
  • 2
    $\begingroup$ have you done a search on the site? there are quite a few questions (and answers) dealing with pure-Neumann conditions for Poisson systems, such as this one or this one. The latter q is concerned with biCGstab for a pure-Neumann problem, which seems very close to what you want. $\endgroup$
    – GoHokies
    Jul 9, 2017 at 12:45
  • 1
    $\begingroup$ You can always post-process your solution to make the constant part vanish; simply subtract the mean. $\endgroup$
    – cfh
    Jul 10, 2017 at 8:49

2 Answers 2

3
$\begingroup$

Krylov solvers have no problem converging if the system is consistent, i.e., if the right-hand side is in the image of the system matrix $A$; see, e.g., Iterative Krylov Methods for Large Linear Systems by van der Vorst. PETSc's CG solver should work.

You can even use a multigrid preconditioner if you make sure that the coarse solver preserves the null space. PETSc's default ILU solver is not sufficient, better use Jacobi; see http://lists.mcs.anl.gov/pipermail/petsc-users/2012-February/012139.html.

$\endgroup$
2
$\begingroup$

The problem I imagine you are trying to solve is a diffusion equation with source term with homogeneous Neumann BC. To do so it must be well posed in order to obtain physical and good results.

The resultant system leads to a singular differential operator matrix $A$, once the BC have been applied. This is so, as you already mentioned, because the following problem (for example I use a linear equation because for the nonlinear case the same would apply ):

$$\left\{\begin{array}{ll}% -\partial_x^2u=f & x\in(0,1)\\ \partial_xu=0 & x=0,x=1 \end{array}\right. \tag{*}$$

is invariant under the transformation $u\to u+constant$.

And therefore some value for $u$ must be specify $\textbf{inside}$ the domain. Many people impose a Dirichlet BC instead one of the Neuman BC causing the solution to differ from the actual one, which now solves:

$$\left\{\begin{array}{ll}% -\partial_x^2u=f & x\in(0,1)\\ u=0 & x=0\\ \partial_xu=0 & x=1\\ \end{array}\right.$$

The problem $(*)$ is well posed if $\int_{0}^{1}{f\,dx}=0$, in fact the source function I propose:

$$f=cos(2\pi x)$$ fits well for the well-posedness of $(*)$, only some reference value for $u$ is required to manage the uniqueness, for example we require that $u(1/2)=1$ $\textbf{without any loss of generality}$. I put these words in bold due to your requirement that your constant must be set to $0$.

N = 100; % #nodes
uref = 1; %Ref value for u for centre node
deltax = 1/(N-1); % step
x = (0:deltax:1)';
f = cos(x*2*pi); %distributed source


b = zeros(N,1); % Source vector
A=zeros(N,N); % Stiffness matrix

for i = 2 : N-1
    A(i,i-1:i+1) = -1/deltax^2*[1, -2, 1];
    b(i) = f(i);
end

%Neumann BC
A(1,1:2) = 2/deltax^2*[1,-1];
b(1) = f(1);
A(N,N-1:N) = 2/deltax^2*[-1, 1];
b(N) = f(N);

%Uniqueness condition e.g. central node
idx = floor(N/2);
A(idx,:) = 0;
A(idx,idx) = 1;
b(idx) = uref;

% System solution
u = A\b;

% Verify the diff equation
ddu = zeros(N,1);
for i = 2 : N-1
     ddu(i) = -(u(i+1)-2*u(i)+u(i-1))/(deltax^2);
end
ddu(1) = f(1);
ddu(N) = f(N);

%Solution residual
res = norm(ddu-f);

%Plots
plot(x,u,'r','linewidth',2); %Solution
figure
plot(x,ddu-f,'k','linewidth',1) % Error

The above code produces for $(*)$ a solution Solution of $(*)$ and the pointwise residual given by $e=\partial_x^2 u+f$ is given below for reference: Residual for $(*)$

You now are free to choose the arbitrary constant (forget the fact that $u(1/2) = 1$) to your problem, just add it.

Formally the uniqueness condition that you mention is simply imposed by the restriction: $$\int_{0}^{1}{u\,dx}=0$$ which constrains $u$ in a way that the transformation $u\to u+constant$ is not valid any more, and therefore the constant must be zero.

This condition would be imposed in an analogous manner when solving the system:

%Uniqueness condition 
idx = N-1;
A(idx,:) = 1;
b(idx) = 0;

Or at the end, when $u_{calc}$ has been obtained, you do: $$u=u_{calc}-\overline{u}_{calc}=u_{calc}-\int_{0}^{1}{u_{calc}\,dx}$$

$\endgroup$
3
  • $\begingroup$ (remark) It is not pointwise error, it is pointwise residual. $\endgroup$
    – VorKir
    Jul 10, 2017 at 17:41
  • $\begingroup$ I thought for a while about your idea - it will work badly in multidimensional case, won't it? You will make a solution constant on a line where actually it is not necessarily constant. $\endgroup$
    – VorKir
    Jul 10, 2017 at 19:13
  • $\begingroup$ No. You only set the value of one point... not of a line. @VorKir $\endgroup$
    – HBR
    Jul 10, 2017 at 21:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.