matlab中消元法方程组,列主消元法解非奇异线性方程组的MATLAB程序

function [x] = LMain_elimination(A,b,prec,n)

%%

% 列主消元法求解线性代数方程组 Ax = b的MATLAB实现

% A为待求解方程组的系数矩阵(要求A为非奇异系数矩阵)

% b为方程组右端常数项(列向量)

% prec为计算过程中的控制常数，当某个选主元或完成消元后的系数小于prec时就停止计算

% n为计算过程中的有效数字长度(用vpa函数保留)

%%

% x为求解的结果

%%

fprintf('列主消元法求解线性方程组 Ax = b');

[col,rol] = size(A); % 求方程组系数矩阵的行数和列数,这里col == rol

M = zeros(col,rol); % 用0初始化M矩阵

x = zeros(col,1); % 用0初始化x向量

%%

% 消元计算:对k = 1,2,…,n-1依次计算

for k = 1:col-1

[maxv,r] = max(abs(A(k:n,k))); % 求第k列的最大值，用maxv保存,r记录的是从第k行起的第r行(最大值所在的行号)，也就是第r + k - 1行

if (maxv < prec)

disp('主元过小，停止计算');

return;

end

r = r + k - 1; % 最大值所在的行号

A([r k],:) = A([k r],:); % 交换系数矩阵的第k行和第r行

b([r k]) = b([k r]); % 交换常数项的第k个数和第r个数(列向量)

if (vpa(A(k,k),n) == 0) % 如果出现除数为0的情况就不必再计算，跳出此函数

return;

else

M(k+1:col,k) = vpa(vpa(A(k+1:col,k),n)/vpa(A(k,k),n),n);

A(k+1:col,k+1:col) = vpa(vpa(A(k+1:col,k+1:col),n) - vpa(M(k+1:col,k),n)*vpa(A(k,k+1:col),n),n);

b(k+1:col) = vpa(vpa(b(k+1:col),n)-vpa(M(k+1:col,k),n)*vpa(b(k),n),n);

end

end

%%

% 回代求解

if (vpa(A(col,rol),n)~=0) % 如果除数不为0,就计算最后一个未知数的根

x(col) = vpa(vpa(b(col),n)/vpa(A(col,rol),n),n);

end

for k = col-1:-1:1

if (vpa(A(k,k),n) == 0) % 除数为0就跳出函数

return;

else

% 依次回代求解其余的根

x(k) = vpa(vpa(vpa(b(k),n) - vpa(A(k,k+1:col),n)*vpa(x(k+1:col),n),n)/vpa(A(k,k),n),n);

end

end

end

%%