1
$\begingroup$

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
$\endgroup$
2
  • 2
    $\begingroup$ 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? $\endgroup$ May 18 at 16:15
  • $\begingroup$ @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? $\endgroup$
    – User124356
    May 18 at 16:50

1 Answer 1

2
$\begingroup$

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.

enter image description here

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.

$\endgroup$
10
  • $\begingroup$ Thank you for your answer. Is there any way to find the solution of ill-posedness ? $\endgroup$
    – User124356
    May 19 at 4:18
  • 1
    $\begingroup$ Let's just think about it in terms of physics: Is there any way to find the solution to the heat conduction problem proposed above? If the source term averages out to zero then there is a solution, otherwise no. You can call it the consequence of the Gauss theorem. $\endgroup$ May 19 at 4:30
  • $\begingroup$ 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$? $\endgroup$
    – User124356
    May 23 at 14:18
  • $\begingroup$ 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. $\endgroup$ May 23 at 14:21
  • $\begingroup$ Do you mean by $u(x_0,y)= u_0$? $\endgroup$
    – User124356
    May 23 at 14:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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