Why are the solutions zeros for a linear equations ?

Hi everyone,

I have a Ax=b need to be solved. A is a 119309*119309 matrix, and b=zeros(119309,1). I used Jacobi and Gauss-seidel methods to solve the equations and the initial values for x are zeros. But all of the solutions for x are zeros, that is zeros(119309,1). I calculate the determinant of A and is 0, that means A is a singular matrix, and the sprank (A)<n(119309). I don't know why I can not get other non-zero solutions. I have been confused for this problem for a week and cannot find the reason. I will appreciate your help if you can help me. My codes in matlab are as follows:

Jacob:

    function[x,k,index]=Jacobi(A,b,ep,it_max)
    step=0
    if nargin<4 
       it_max=100;
    end
    if nargin<3 
       ep=1e-10;
    end
    n=length(A);
    k=0;
    x=zeros(n,1);
    y=zeros(n,1);
    index=1;
    while 1
        for i=1:n
            y(i)=b(i);
            for j=1:n
                if j~=i 
                    y(i)=y(i)-A(i,j)*x(j);
                end
            end
            if abs(A(i,i))<1e-10 | k==it_max
                index=0;
                return;
            end
            y(i)=y(i)/A(i,i);
        end
        if norm(y-x,inf)<ep
            break;
        end
        x=y;
        k=k+1;
        step=step+1

    end

    [x,k,index]=Jacobi(A,b,1e-10,10)

Gauss-seidel:

    function[v,sN,vChain]=GaussSeidell(A,b,x0,errorBound,maxSp)
    step=0
    error=inf;
    vChain=zeros(15,151);
    k=1;
    fx0=x0;
    L=-tril(A,-1);
    U=-triu(A,1);
    C=diag(diag(A));
    b=zeros(151,1);
    while error>=errorBound & step<maxSp
        x0=inv(C)*(L+U)*x0+inv(C)*b;
        vChain(k,:)=x0;
        k=k+1;
        error=norm(x0-fx0);
        fx0=x0;
        step=step+1;
        e=1
    end
    v=x0;
    sN=step

    end

    [v,sN,vChain]=GaussSeidell(A,b,x0,1e-6,100)

[Pierre Benoit]

First, please format your code.

One way to see why a program's behaviour is wrong is to test a small example and see what it does step by step (with inserting breakpoints for example). You could start with a matrix A 2x2 like

A = [1 2; 1 2];

where one of the solution is

x = [-2 ; 1]

And from what I gathered of the methods you're using to solve this problem, you use a starting point and iterates, but here you start with an obvious solution which is 0, so maybe your code just stop there since you already found a solution. You may want to start with a different vector x.

Answer 1

[Pierre Benoit]

The null function in Matlab achieve what you desire (for b=0) but I suppose this is homework related to better understand these algorithms rather than an efficient implementation to solve this.

As I stated in my comment, if you start with a solution, your algorithm will end in the first iteration, so you have to start by another vector. And you have to remember the fact that such an algorithm can never converge or converge quite slowly. Trials and errors may be necessary to find a good starting point.

Also, are you sure you can use these methods to solve Ax = 0 ? In dimension 2, you will always come back to the same vector every two iteration since det(A) = 0. I don't know for higher dimensions but all the information I found was about to solve Ax=b with b different from 0.

Remember that such equation (b=0) have infinite solutions for A singular, see Kernel for more informations.

[rong]

Yes, I changed the initial solution x=0.0001 and tried for a smaller matrix equations. Then I got a different solution with 0. But when I applied on the original matrix A(119309*119309). Errors just displayed: Warning: "Matrix is singular to working precision." or Error"using inv Out of memory. Type HELP MEMORY for your options. "

Answer 2

[Matt J]

In every pass through the loop, you reset the state variables to some unchanging quantity. In the Jacobi, you have

    y(i)=b(i);

while in the Gauss Seidel, you have

    x0=inv(C)*(L+U)*x0+inv(C)*b;

The algorithm can't progress if you're always resetting...

Answer 3

[rong]

Does Jacobi or Gauss-seidel solve Ax=b if A is singular matrix? If not I tried x=pinv(A)*b, but it only works for full matrix. In my case, A is a sparse matrix, so I used QR to decompose A, but I am not sure if it is correct?

[Matt J]

Why not just do A\b? For example,

    >> A=sparse([ones(3,2),zeros(3,5)]);  b=[1;1;1];
    >> x=A\b
    Warning: Rank deficient, rank = 1, tol =  7.691851e-14. 

    x =

         1
         0
         0
         0
         0
         0
         0