A question on the Poisson equation with Neumann and periodic boundary conditions on a rectangular region

Question

I am trying to solve the following PDE by using finite difference

\begin{eqnarray*} \Delta u&=& f ~~on~~(0,1)\times(0,1)\\ \frac{\partial u}{\partial y}(x,0)&=&0=\frac{\partial u}{\partial y}(x,1)~for ~x\in[0,1]\\ u(0,y)&=&u(1,y)~for ~y\in[0,1]\\ \end{eqnarray*}

For a uniform spacing $h$, I got the following equation, $$ \frac{1}{h^2}(u_{{i-1},j}-4u_{i,j}+u_{{i+1},j}+u_{{i},{j-1}}+u_{{i},{j+1}})=f_{i,j} $$ for $i =1,2,......,Nx+2$ and $j=1,2,...,Ny+2$. After the implementation of the boundary conditions, I converted the system into the system of linear equations $$ Au=f $$.

Now, I am trying to solve this system of linear equations, for certain value of $Nx$, I got the following warning;

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.053110e-20.

I believe this message is appearing because the matrix $A$ become singular after the implementation of the boundary conditions. The following a sub-code of my main code to solve this problem.

function [v1] = new_v1(w,Nx,Ny,dx,dy)
% -------------------------------------------------------
Iint = 1:Nx+1; % Indices of interior point in x direction
Jint = 1:Ny+2; % Indices of interior point in y direction
%---------------------------------------------
%assembly of the tridiagonal matrix(LHS)
sx = 1/(dx^2);
sy = 1/(dy^2);
e=ones(Nx+1,1);
T=spdiags([sx.e,((-2sx)+(-2sy)).e,sx.*e],-1:1,Nx+1,Nx+1);
T(1,Nx+1)= sx;
T(Nx+1,1)= sx;
D=spdiags(sy.*e,0,Nx+1,Nx+1);
A=blktridiag(T,D,D,Ny+2);
for i=1:Nx+1
    for j=1:Nx+1
        A(i,j)=(1/2).A(i,j);
        A((Nx+1)(Ny+1)+i,(Nx+1)(Ny+1)+j)=(1/2).A((Nx+1)(Ny+1)+i,(Nx+1)(Ny+1)+j);
    end
end
%---------------------------------------------------------------
%Solve the linear system
rhs = w ;
for i=1:Nx+1
rhs(i,1)=(1/2).*rhs(i,1);
rhs(i,Ny+2)=(1/2).*rhs(i,Ny+2);
end

%convert the rhs into column vector
F = reshape(rhs, (Nx+1)*(Ny+2),1);

uvec = A\F;
v1(Iint, Jint)= reshape(uvec, Nx+1,Ny+2);

end

Answer 1

Let's give this PDE a simple physical interpretation as the problem of finding a steady-state temperature distribution on the side surface of a cylinder.

The x coordinate describes the length along the azimuthal direction (normalized to unity), and the y coordinate describes the length along the cylinder axis (normalized to unity). The variable $u(x,y)$ is the temperature, and the right-hand side $f$ is the heat source, and the heat diffusivity coefficient is equal to 1. The condition of $\partial_y u = 0$ at the top and bottom of the cylinder is satisfied because those are the boundaries of the domain, so the heat flux normal to the boundary vanishes there. The condition $u(0,y)=u(1,y)$ is naturally satisfied because the system is periodic in $x$. This physics problem in general is not solvable, there is no steady-state solution here unless the integral of the heat source over the surface vanishes, $\int \! \!{f}dxdy$=0. In the latter case, there is a smooth solution defined up to a constant.

The BC $u(0,y)=u(1,y)$ is less stringent than the periodicity conditions in the physical system where $u_x(0,y)=u_x(1,y)$ is also imposed (assuming $f$ does not contain delta-functions). To account for that, one can put the surface of the cylinder in contact with a thermostat along the line $x=0$, thus enforcing an arbitrary temperature $T_0(y)$ along this line, and that would still be a solution satisfying the PDE and all given boundary conditions. Therefore there is an infinite number of solutions to this problem, corresponding to the choice of $T_0(y)$, so the problem as stated is ill posed. This ill-posedness is manifested in the matrix $A$ singularity. By imposing additional constraints on the problem one could make it well posed.