M010K
0%

计算方法实验/实验一

发表于 更新于 分类于 阅读次数： Valine：
本文字数： 2.5k 阅读时长 ≈ 2 分钟
解线性方程组的高斯消去法与高斯列主元消去法

实验一

解线性方程组的高斯消去法与高斯列主元消去法

【作业内容】：解线性方程组：

  1. 高斯消去法

    用高斯消去法求解上述方程组，输出消元过程以及$Ax=b$中矩阵$A$、$b$、解向量$x$及$detA$。

  2. 高斯列主元消去法

    用高斯列主元消去法求解上述方程组，输出列主元选取、消元过程，包括行交换次序，$Ax=b$中矩阵$A$、$b$、解向量$x$及$detA$。

  3. 比较两种消去法. 用MATLAB的内部函数$inv$求出系数矩阵的逆矩阵，再输入命令$x=inv(A)*b$, 即可求出上述两个方程组的解。并与高斯消去法、高斯列主元消去法求出的解进行比较，体会选主元的方法具有良好的数值稳定性。

说明：计算过程保留小数点以后8位数字

高斯消去法

  • matlab代码如下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
% Gauss.m
% 高斯法求解线性方程组Ax=b
% A为输入矩阵系数,b为方程组右端系数
% 方程组的解保存在x变量中
% 先输入方程系数

% 设置计算过程的默认精度
digits(9);

A = vpa([3.01,6.03,1.96; 1.348,6.65,0.002; 1.47, -4.81,9.34])
disp('detA=')
disp(det(A))
b = [11,8,6]'
[m,n]=size(A);
%检查系数正确性
if m~=n
    error('矩阵A的行数和列数必须相同');
    return;
end
if m~=size(b)
    error('b的大小必须和A的行数或A的列数相同');
    return;
end
%再检查方程是否存在唯一解
if rank(A)~=rank([A,b])
    error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');
    return;
end

%这里采用增广矩阵行变换的方式求解
c=n+1;
A(:,c)=b; 


%高斯消元，输出消元过程
disp('高斯消元:');
for k=1:n-1
    A(k+1:n, k:c) = A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c)
end

%%回代结果
x=zeros(length(b),1);  
x(n)=A(n,c)/A(n,n);
for k=n-1:-1:1
x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

%显示计算结果
disp('x=');
disp(x);

高斯列主元消去法

  • matlab代码如下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
% Column_Gauss.m
% 高斯法求解线性方程组Ax=b
% A为输入矩阵系数,b为方程组右端系数
% 方程组的解保存在x变量中
% 先输入方程系数

% 设置计算过程的默认精度
digits(9);
A = vpa([3.01,6.03,1.96; 1.348,6.65,0.002; 1.47, -4.81,9.34])
disp('detA=')
disp(det(A))
b = [11,8,6]'
[m,n]=size(A);

%检查系数正确性
if m~=n
    error('矩阵A的行数和列数必须相同');
    return;
end
if m~=size(b)
    error('b的大小必须和A的行数或A的列数相同');
    return;
end
%再检查方程是否存在唯一解
if rank(A)~=rank([A,b])
    error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');
    return;
end

%这里采用增广矩阵行变换的方式求解
c=n+1;
A(:,c)=b; 

for i=1:n
    disp(['第',num2str(i),'次交换列主元前:'])
    disp(A);
    [max_col,position_to_i]= max(abs(A(i:end,i)));%只考虑i后的行列，选取绝对值较大的作为主元
    if max_col == 0
       error('主对角元是0，无唯一解，无法用高斯主列消元法解');
    end
    %交换最大元列 的行
    max_col_rowposition = position_to_i+i-1;%列中的实际位置..
    if max_col_rowposition ~= i
        temp = A(i,:);
        A(i,:) = A(max_col_rowposition,:);
        A(max_col_rowposition,:) = temp;
    end
    disp(['第',num2str(i),'次交换列主元后:'])
    disp(A);
    for j = (i+1):n
        A(j,:) = A(j,:) - A(j,i) / A(i,i) * A(i,:);
    end
    disp(['第',num2str(i),'次消元后:'])
    disp(A);
    
end


%%回代结果
x=zeros(length(b),1);  
x(n)=A(n,c)/A(n,n);
for k=n-1:-1:1
    x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

%显示计算结果
disp('x=');
disp(x);

Matlab逆矩阵求解

1
2
3
4
A = [3.01,6.03,1.96; 1.348,6.65,0.002; 1.47, -4.81,9.34];
b = [11,8,6]';
x = inv(A)*b;
disp(x);