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

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

For a uniform spacing $h$, I got the following equation,

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

[Sharmin Kibria]

Hi kaps,

I believe the warning is coming from the matrix inversion uvec = A\F;. These warnings are generally caused when the singular value of the matrix A becomes very close to zero. In that case the determinant of the matrix becomes very close to zero.

Sometimes these errors are fixed by modifying the singular values to ensure they do not becomes close to numerically zero.

[U,S,V] = svd(A) performs a singular value decomposition of matrix A, such that A = U*S*V'. Setting the singular values to (S+e) (e is a very small number) instead of S should resolve the issue.

Thanks

Sharmin

Answer 1

[Christine Tobler]

Yes, the matrix A becomes singular after applying the last two for-loops (which I think are the boundary conditions).

You can verify that by calling null(full(A)) on an example matrix (I used Nx = Ny = 10, dx = dy = 0.1). This showed that there is a null space of dimension one, and the vector in that null space had all elements of equal value.

That means that if you insert all values of x as being the same number, then A*x = 0 is always true. That means there is a degree of freedom left in this linear system.

Thinking about your initial equations, this makes sense: None of the three equations above tie the value of u to any specific number: The derivatives are specified, but we just know that u(0, y) = u(1, y) is required, which could be true for any function u that is shifted by a constant.

I would suggest just setting one of the values of x to an arbitrary value, and then rearranging the rest of this linear system based on that value.

[kaps]

Do you mean I need to choose some specific value at the boundary, so I get unique solution? I know that if u is the solution of the given problem then u+c is also a solution of the given problem.

[Christine Tobler]

Yes, when there is no unique solution, the linear system is singular, a warning is given, and the results are less reliable.

You can make the solution unique by adding an additional equation which sets any one of the values in the solution to zero (or any other value, basically this fixes the c). To keep the linear system square (which is cheaper to solve typically), you could remove one of the existing equations, since the system being singular means that one of these equations is redundant.

Or you could use the lsqminnorm function, which computes the solution x which has minimal 2-norm among all solutions. But keep in mind that this can get expensive as the size of the system increases.