Geek Answers Handbook

Gauss exception in Matlab

I use the matlab code from this book: http://books.google.com/books/about/Probability_Markov_chains_queues_and_sim.html?id=HdAQdzAjl60C Here is the code:

    function [pi] = GE(Q)

    A = Q';
    n = size(A);
    for i=1:n-1
      for j=i+1:n
         A(j,i) = -A(j,i)/A(i,i);
      end
         for j =i+1:n
            for k=i+1:n
        A(j,k) = A(j,k)+ A(j,i) * A(i,k);
         end
      end
      end

     x(n) = 1;
      for i = n-1:-1:1
        for j= i+1:n
          x(i) = x(i) + A(i,j)*x(j);
        end
       x(i) = -x(i)/A(i,i);
      end

      pi = x/norm(x,1);

Is there a faster code that I don't know about? I call this function millions of times, and it takes too much time.

+5
source share
6 answers

MATLAB has a whole set of built-in linear algebra routines - type help slash, help luor help chol, to start with several common ways to efficiently solve linear equations in MATLAB.

LAPACK/BLAS, , , . "" , MATLAB, , , m .

, .

+9

, Matlab (mldivide) , , , lu. , mldivide , (, ).

, mldivide lu, C Fortran, Matlab . , , , , (, ).

: pivoting, , , , .

, O (n 3). - , -, .

+8
function x = naiv_gauss(A,b);
n = length(b); x = zeros(n,1);
for k=1:n-1 % forward elimination
      for i=k+1:n
           xmult = A(i,k)/A(k,k);
           for j=k+1:n
             A(i,j) = A(i,j)-xmult*A(k,j);
           end
           b(i) = b(i)-xmult*b(k);
      end
end
 % back substitution
x(n) = b(n)/A(n,n);
for i=n-1:-1:1
   sum = b(i);
   for j=i+1:n
     sum = sum-A(i,j)*x(j);
   end
   x(i) = sum/A(i,i);
end
end
+5

, Ax = d A d - . "A" "LU", "LU-", matlab, : LUx = d Matlab : [L, U] = lu (A) U L , A = LU. L . (https://www.mathworks.com/help/matlab/ref/lu.html)

, , Ly = d, y = Ux. x , , , y, x : = \; = U\

x.

, (.. A d ), , , .

, , .

+1

(aka ) n n :

function A = naiveGauss(A)

% find the size

n = size(A);
n = n(1);

B = zeros(n,1);

% We have 3 steps for a 4x4 matrix so we have
% n-1 steps for an nxn matrix
for k = 1 : n-1     
    for i = k+1 : n
        % step 1: Create multiples that would make the top left 1
        % printf("multi = %d / %d\n", A(i,k), A(k,k), A(i,k)/A(k,k) )
        for j = k : n
            A(i,j) = A(i,j) - (A(i,k)/A(k,k)) *  A(k,j);
        end
        B(i) = B(i) - (A(i,k)/A(k,k))  * B(k);
    end
end
0
function Sol = GaussianElimination(A,b)


[i,j] = size(A);

for j = 1:i-1


    for i = j+1:i


        Sol(i,j) = Sol(i,:) -( Sol(i,j)/(Sol(j,j)*Sol(j,:)));


    end


end


disp(Sol);


end
0

More articles:


All Articles