eliminación de Gauss en Matlab

Question

Estoy utilizando el código Matlab de este libro: http://books.google.com/books/about/Probability_Markov_chains_queues_and_sim.html?id=HdAQdzAjl60C Aquí está el código:

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); 

¿Hay un código más rápido que yo no conozco? Llamo a esto funciones millones de veces y lleva demasiado tiempo.

Answer 1

MATLAB tiene un conjunto de rutinas incorporadas de álgebra lineal - Tipo help slash, help lu o help chol a empezar con algunas de las formas más comunes para resolver ecuaciones lineales de manera eficiente en MATLAB.

Debajo del capó, estas funciones generalmente llaman a las rutinas de biblioteca LAPACK/BLAS optimizadas, que generalmente son la forma más rápida de hacer álgebra lineal en cualquier lenguaje de programación. Comparado con un lenguaje "lento" como MATLAB, no sería inesperado si fueran órdenes de magnitud más rápidos que una implementación de m-file.

Espero que esto ayude.

Answer 2

A menos que específicamente desee implementar su propia, debe utilizar el operador de barra invertida de Matlab (mldivide) o, si desea los factores, lu. Tenga en cuenta que mldivide puede hacer más que la eliminación gaussiana (por ejemplo, hace mínimos cuadrados lineales, cuando corresponda).

Los algoritmos utilizados por mldivide y lu son de las bibliotecas C y Fortran, y su implementación en Matlab nunca será tan rápida. Sin embargo, si está decidido a usar su propia implementación y desea que sea más rápida, una opción es buscar formas de vectorizar su implementación (quizás inicie here).

Otra cosa a tener en cuenta: la implementación de la pregunta no hace ningún pivoting, por lo que su estabilidad numérica generalmente será peor que una implementación que pivotea, e incluso fallará en algunas matrices no singulares. existen

diferentes variantes de eliminación de Gauss, pero todos ellos son O (n) algoritmos. Si un enfoque es mejor que otro depende de su situación particular y es algo que necesitaría investigar más.

Answer 3

Para el enfoque ingenuo (aka sin fila swapping) para un n por n matriz:

function A = naiveGauss(A) 

% find's 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 

Answer 4

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 

Answer 5

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 

Answer 6

Asumamos Ax = d donde a y d se conocen matrices. queremos representar "A" como "L U" utilizando la función de "descomposición LU" incrustado en Matlab así: LUx = d Esto se puede hacer en Matlab siguiente: [L, U] = lu (A) que en términos devuelve una matriz triangular superior en U y una matriz triangular inferior permutada en L tal que A = L U. El valor de retorno L es un producto de matrices triangulares y de permutación inferiores. (https://www.mathworks.com/help/matlab/ref/lu.html)

Entonces, si asumimos Ly = d donde y = Ux. Como x es Desconocido, por lo tanto, y es desconocido también, al conocer y encontramos x de la siguiente manera: y = L \ d; x = U \ y

y la solución se almacena en x.

Esta es la forma más simple de resolver el sistema de ecuaciones lineales siempre que las matrices no sean singulares (es decir, el determinante de la matriz A yd no sea cero); de lo contrario, la calidad de la solución no sería tan buena como se esperaba y puede dar resultados incorrectos.

si las matrices son singulares, por lo tanto, no pueden invertirse, se debe utilizar otro método para resolver el sistema de las ecuaciones lineales.