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

17 views (last 30 days)
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
  1 Comment
Sharmin Kibria
Sharmin Kibria on 23 May 2022
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

Sign in to comment.

Answers (1)

Christine Tobler
Christine Tobler on 25 May 2022
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.
  2 Comments
Christine Tobler
Christine Tobler on 27 May 2022
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by