高斯消去法,通常包含两个过程,消去过程和回代过程,不同的高斯消去法其消去过程也不尽相同,一般的高斯消去法旨在将系数阵通过行变换,变成上三角矩阵,而列主元素的高斯消去法,在一般的基础上,每次行变换时将该列绝对值最大元素变换到主元位置上,而量化的列主元素高斯消去法则不进行行交换,通过一个l向量来存储列主元素所在行的位置,按照该顺序进行行变换,可以减少用于行交换的处理时间。
直接上代码,matlab编写:
1.高斯消去法
function x = gauss1(A,b) %A为系数阵,b为常数向量
x1 = A;
x2 = b;
size_x1 = size(x1);
n = size_x1(1); %n为方程数
if A(1,1) == 0
x = 'error';
disp('error!!!');
return;
end %对于对角元为0不适用
for i = 1:1:n
col_i = x1(:,i); %取列向量
if col_i(i) == 0
x = 'error';
disp('error!!!');
return;
end
m_i = col_i / col_i(i);
for k1 = i+1:1:n
x2(k1) = x2(k1) - m_i(k1) * x2(i);
for k2 = 1:1:n
x1(k1,k2) = x1(k1,k2) - m_i(k1) * x1(i,k2);
end
end
end
%回代过程
res_x1 = x1;
res_x2 = x2;
x = zeros(size(res_x2));
for i = n:-1:1
sum = 0;
if i == n
x(n) = x2(n) / x1(n,n);
else
temp = res_x1(i,:) .* x';
for j = 1:1:n
sum = sum + temp(j);
end
x(i) = (x2(i) - sum) / x1(i,i);
end
end
return;
end
2.列主元素高斯消去法
function x = gauss2(A,b)
x1 = A;
x2 = b;
size_x1 = size(x1);
n = size_x1(1);
for i = 1:1:n
col_i = x1(:,i);
for ii = i+1:1:n
if col_i(ii) > col_i(i)
t1 = zeros(size(col_i));
for jj = 1:1:n
t1(jj) = x1(i,jj);
x1(i,jj) = x1(ii,jj);
x1(ii,jj) = t1(jj);
end
t2 = x2(i);
x2(i) = x2(ii);
x2(ii) = t2;
end
col_i = x1(:,i);
end
if col_i(i) == 0
x = 'error';
disp('error!!!');
return;
end
m_i = col_i / col_i(i);
for k1 = i+1:1:n
x2(k1) = x2(k1) - m_i(k1) * x2(i);
for k2 = 1:1:n
x1(k1,k2) = x1(k1,k2) - m_i(k1) * x1(i,k2);
end
end
end
res_x1 = x1;
res_x2 = x2;
x = zeros(size(res_x2));
for i = n:-1:1
sum = 0;
if i == n
x(n) = x2(n) / x1(n,n);
else
temp = res_x1(i,:) .* x';
for j = 1:1:n
sum = sum + temp(j);
end
x(i) = (x2(i) - sum) / x1(i,i);
end
end
return;
end
3.量化的列主元素高斯消去法
function x = gauss3(A,b)
size_A = size(A);
n = size_A(1);
l = 1:1:n;
S = zeros(size(l));
for i = 1:1:n
for j = 1:1:n
if abs(A(i,j)) > S(i)
S(i) = abs(A(i,j));
end
end
end
for i = 1:1:n
max_as = 0;
col_i = A(:,i);
for j = i:1:n
if col_i(l(j)) / S(j) > max_as
max_as = col_i(l(j)) / S(j);
t = l(i);
l(i) = l(j);
l(j) = t;
end
end
m_i = col_i / col_i(l(i));
for k1 = i+1:1:n
b(l(k1)) = b(l(k1)) - m_i(l(k1)) * b(l(i));
for k2 = 1:1:n
A(l(k1),k2) = A(l(k1),k2) - m_i(l(k1)) * A(l(i),k2);
end
end
end
x = zeros(size(b));
for i = n:-1:1
sum = 0;
if l(i) == l(n)
x(n) = b(l(n)) / A(l(n),n);
else
temp = A(l(i),:) .* x';
for j = 1:1:n
sum = sum + temp(j);
end
x(i) = (b(l(i)) - sum) / A(l(i),i);
end
end
return;
end