gauss.m

数值分析上机实验

function[x,det,index] = gauss(a,b)
[n,m] = size(a);
nb = length(b);
if n ~= m
  error('the row and columns of matrix a must be equal!');
return;
end
if m ~= nb
  error('the columns of a must be length of b');
return;
end

index = 1;
det = 1;
x = zeros(n,1);
for k = 1 : (n-1)
  a_max = 0;
for i = k:n
  if abs(a(i,k)) > a_max
    a_max = abs(a(i,k));
    r = i;
end
end
if a_max < 1e-10
  index = 0;
return;
end

if r>k
  for j = k:n
    z = b(k);
    b(k) = b(r);
    b(r) = z;
det = -det;
end

for i = k+1:n
  m = a(i,k)/a(k,k);
  for j = k+1:n
    a(i,j) = a(i,j) - m*a(k,j);
  end
  det = det*a(k,k);
end
det = det*a(n,n);
if abs(a(n,n)) < 1e - 10
  index = 0;
 return;
end

for k = n:-1:1
  for j = k+1:n
    b(k) = b(k) - a(k,j) * x(j);
end
x(k) = b(k)/a(k,k);
end