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? May 18 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? May 18 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$? May 23 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. May 23 at 14:21
• Do you mean by $u(x_0,y)= u_0$? May 23 at 14:26