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.

[John D'Errico]

As Christine has said, that is exactly the source of your problem. You admit that you can essentially translate the problem by any constant value, and the solution is as good. This means you need to choose any arbitrary point, and force the solution to pass through that point. If you don't, then the problem becomes singular. And backslash won't solve a singular problem. (At least not in a way that will make you happy.)

Also as Christine has suggested, one solution is to use lsqminnorm, which will be expensive. But lsqminnorm actually does implicitly choose some constant for you, It just does so in a way that won't be terribly useful to you. Instead, the correct way to solve the problem is to indeed enforce that constant to take on some specific value. Just pick one corner of the domain, and set u(0,0)=0, for example as another boundary condition. That will immediately resolve the problem.

If you think about it, a singular system of equations in linear algebra typically means there are infinitely many equivalent solutions. The singularity just reflects back to you what you already know. But it is also what gives the solver a case of the hiccups, since it will not resolve that singularity for you.

This is a common issue faced by students in engineering courses. They might be told to build a truss model of a bridge that minimizes the potential energy of the structure. But unless they build in a boundary condition that properly fixes the structure at some location in space, the equations become singular. Essentially the entire structure can move as a whole anywhere in the universe, as long as it has the same fundamental shape.