Sunday, October 11, 2009

实现最简单的Gaussian Elimination

不考虑任何特殊情况,最简单的高斯消去法可实现如下,

% the most simple Gauss Elimination

function b = Gauss_Elim(A,b)

n = length(b);

%Elimination phase
for i =1:n-1
for j =i+1:n
if(A(j,i))
l = A(j,i)/A(i,i);
A(j,i:n)= A(j,i:n)- l*A(i,i:n);
b(j) = b(j)-l*b(i);
end
end
end

%back substitution
b(n)= b(n)/A(n,n);
for i = (n-1):-1:1
b(i)=b(i)-A(i,i+1:n)*b(i+1:n);
b(i)=b(i)/A(i,i);
end

试比较下面的代码,

function b = g(A,b)
n = length(b);
%-----------------Elimination phase-------------
for k = 1:n-1
for i = k+1:n
if A(i,k) ~= 0
lambda = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - lambda*A(k,k+1:n);
b(i)= b(i) - lambda*b(k);
end
end
end

size(A)

%--------------Back substitution phase-----------
for k = n:-1:1
k
b(k) = (b(k) - A(k,k+1:n)*b(k+1:n))/A(k,k)
end

这两段代码的区别,
1 在消去阶段,区别在于对主对角线元素下元素的处理。前者作了处理,后者未作处理。
2在回代阶段,后者代码中当k=n时,A(k,k+1:n)相当于A(n,n+1:n),在matlab中实际上是一个Empty matrix: 1-by-0,
但是当A(k,k+1:n)*b(k+1:n),结果会是0。这可能是matlab内部的一个设计,使得代码变得简洁。

No comments: