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

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,((-2*sx)+(-2*sy)).*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

• Let's forget the discretization entirely for a sec. Suppose you have a putative solution $u$ to this boundary value problem; for any constant $c$, the function $u + c$ is also a solution to your boundary value problem. What are the implications of this fact and what do you think it means for how you discretize the problem? Commented May 18, 2022 at 16:15
• @DanielShapero this differential equation has infinitely many solutions. So the rank of the matrix A is less than the order of A. Could you please give me little more details so I can think about your point here? Commented May 18, 2022 at 16:50

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.
• I understood your argument. I used the condition, $u_x(0,y)=u_x(1,y)$ and took the source term is zero. But my matrix $A$ is again a singular matrix. I used central difference for the condition $u_x(0,y)=u_x(1,y)$. Could you please write your matrix $A$? Commented May 23, 2022 at 14:18
• Even with the derivative matching conditions, the solution is still defined up to a constant. We need to use one more constraint, $u(x_0,y_0)=u_0$, i.e., enforcing $u$ to be some prescribed constant at a given point. That will fix the matrix singularity. Commented May 23, 2022 at 14:21
• Do you mean by $u(x_0,y)= u_0$? Commented May 23, 2022 at 14:26