第三次MATLAB程序实训
方程组：
{ 2 x ! − x 2 + 3 x 3 = 1 4 x 1 + 2 x 2 + 5 x 3 = 4 x 1 + 2 x 2 = 7 \begin{cases} 2x_! - x_2 + 3x_3 = 1\\4x_1 +2x_2 +5x_3 = 4\\ x_1+2x_2=7\\ \end{cases} ⎩ ⎨ ⎧2x!−x2+3x3=14x1+2x2+5x3=4x1+2x2=7
程序运行图
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-em6jfrHO-1664507204070)(C:\Users\86131\AppData\Roaming\Typora\typora-user-images\image-20220930103447808.png)]
代码设计
%程序缺陷：未判断程序输入是否有效，以及输入的矩阵是否能通过高斯消元法求解
clc
A=input('请输入系数矩阵（方阵A）：');
b=input('请输入常数向量（列向量b）：');
H = size([A b]);
[Ax,bx]=GSelimination_NoChange(A,b,H(1));
[x]=GSback(Ax,bx,H);
disp("(未交换主元）Ax=b的解向量x：")
%format long
disp(x);
[Ax,bx]=GSelimination_Change(A,b,H(1));
[x]=GSback(Ax,bx,H);
disp("(交换主元）Ax=b的解向量x：")
%format long
disp(x)
function [Ax,bx] = GSelimination_Change(A,b,n)
disp("交换主元的高斯消元法:")
%消去过程
%n表示系数矩阵(方阵)的规模
%A表示系数矩阵
%b表示常数矩阵
%消去前先交换主元
for k = 1 : n - 1
ip = k;
for i = k + 1 : n % 寻找主元的下标
if abs(A(i,k)) >= abs(A(k,k))
ip = i;
end
end
if(ip ~= k)
A([k ip],:) = A([ip k],:); %交换行
b([k ip],:) = b([ip k],:);
end
disp([A b])
if(A(k,k)~=0) %主元素不能为0,
for i=(k+1):n %循环行
c=(-1*A(i,k))/A(k,k); %倍乘因子
for j=1:n %单行循环
A(i,j)=A(i,j)+c*A(k,j); %更新矩阵元素
end
b(i)=b(i)+c*b(k); %更新常数向量
end
else
disp("主元素出现零，程序错误！"); %不更换主元，高斯消元并不总能进行
return;
end
end
disp([A b])
Ax=A;bx=b;
end
function [Ax,bx]=GSelimination_NoChange(A,b,n)
disp("不交换主元的高斯消元法：")
%消去过程
%n表示系数矩阵(方阵)的规模
%A表示系数矩阵
%b表示常数矩阵
for k=1:(n-1) %要进行n-1次消去过程
disp([A b])
if(A(k,k)~=0) %主元素不能为0,
for i=(k+1):n %循环行
c=(-1*A(i,k))/A(k,k); %倍乘因子
for j=1:n %单行循环
A(i,j)=A(i,j)+c*A(k,j); %更新矩阵元素
end
b(i)=b(i)+c*b(k); %更新常数向量
end
else
disp("主元素出现零，程序错误！"); %不更换主元，高斯消元并不总能进行
return;
end
end
disp([A b])
Ax=A;bx=b;
end
function [x]=GSback(Ax,bx,n)
%回代过程
%Ax为经过消去过程后得到的上三角矩阵
%bx为经过消去过程后得到的常数向量
%n表示系数矩阵的规模
for i=n:-1:1
if(Ax(i,i)~=0)
x(i)=bx(i);
for j=n:-1:i+1
x(i)=x(i)-Ax(i,j)*x(j);
end
x(i)=x(i)/Ax(i,i);
else
disp("主元素出现零，程序错误！");
return;
end
end
x=x';
end
运行结果
请输入系数矩阵（方阵A）：[2,-1,3;4,2,5;1,2,0]
请输入常数向量（列向量b）：[1;4;7]
不交换主元的高斯消元法：
2 -1 3 1
4 2 5 4
1 2 0 7
2.000000000000000 -1.000000000000000 3.000000000000000 1.000000000000000
0 4.000000000000000 -1.000000000000000 2.000000000000000
0 2.500000000000000 -1.500000000000000 6.500000000000000
2.000000000000000 -1.000000000000000 3.000000000000000 1.000000000000000
0 4.000000000000000 -1.000000000000000 2.000000000000000
0 0 -0.875000000000000 5.250000000000000
(未交换主元）Ax=b的解向量x：
9
-1
-6
交换主元的高斯消元法:
4 2 5 4
2 -1 3 1
1 2 0 7
4.000000000000000 2.000000000000000 5.000000000000000 4.000000000000000
0 -2.000000000000000 0.500000000000000 -1.000000000000000
0 1.500000000000000 -1.250000000000000 6.000000000000000
4.000000000000000 2.000000000000000 5.000000000000000 4.000000000000000
0 -2.000000000000000 0.500000000000000 -1.000000000000000
0 0 -0.875000000000000 5.250000000000000
(交换主元）Ax=b的解向量x：
9
-1
-6
>>
由此可知，对于本例中过程中主元虽不同，但没有出现主元接近于0的情况，所以交换和不交换主元对于结果没有太大影响，两种方法所求结果相同。
方程组二：
{ 8.77 B + 2.40 C + 5.66 D + 1.55 E + 1.0 F = − 32.04 4.93 B + 1.21 C + 4.48 D + 1.10 E + 1.0 F = − 20.07 3.53 B + 1.46 C + 2.92 D + 1.21 E + . 10 F = − 8.53 5.05 B + 4.04 C + 2.51 D + 2.01 E + 1.0 F = − 1.04 3.54 B + 1.04 C + 3.47 D + 1.02 E + 1.0 F = − 12.04 \begin{cases} 8.77B + 2.40C + 5.66D + 1.55E + 1.0F = -32.04\\4.93B + 1.21 C+4.48D+1.10E+1.0F=-20.07\\ 3.53B +1.46C+2.92D+1.21E+.10F=-8.53\\5.05B+4.04C+2.51D+2.01E+1.0F = -1.04\\3.54B+1.04C+3.47D+1.02E+1.0F=-12.04\\ \end{cases} ⎩ ⎨ ⎧8.77B+2.40C+5.66D+1.55E+1.0F=−32.044.93B+1.21C+4.48D+1.10E+1.0F=−20.073.53B+1.46C+2.92D+1.21E+.10F=−8.535.05B+4.04C+2.51D+2.01E+1.0F=−1.043.54B+1.04C+3.47D+1.02E+1.0F=−12.04
运行结果
请输入系数矩阵（方阵A）：[8.77 2.40 5.66 1.55 1.0;4.93 1.21 4.48 1.10 1.0;3.53 1.46 2.92 1.21 1.0;5.05 4.04 2.51 2.01 1.0;3.54 1.04 3.47 1.02 1.0]
请输入常数向量（列向量b）：[-32.04;-20.07;-8.53;-6.30;-12.04]
不交换主元的高斯消元法：
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
4.930000000000000 1.210000000000000 4.480000000000000 1.100000000000000 1.000000000000000
3.530000000000000 1.460000000000000 2.920000000000000 1.210000000000000 1.000000000000000
5.050000000000000 4.040000000000000 2.510000000000000 2.010000000000000 1.000000000000000
3.540000000000000 1.040000000000000 3.470000000000000 1.020000000000000 1.000000000000000
6 列
-32.039999999999999
-20.070000000000000
-8.529999999999999
-6.300000000000000
-12.039999999999999
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 -0.139144811858609 1.298266818700114 0.228677309007982 0.437856328392246
-0.000000000000000 0.493979475484607 0.641801596351197 0.586111744583808 0.597491448118586
0 2.658015963511973 -0.749179019384265 1.117468643101482 0.424173318129989
0 0.071242873432155 1.185347776510832 0.394344355758267 0.596351197263398
6 列
-32.039999999999999
-2.058916761687573
4.366374002280503
12.149486887115163
0.892907639680731
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 -0.139144811858609 1.298266818700114 0.228677309007982 0.437856328392246
-0.000000000000000 -0.000000000000000 5.250792428091452 1.397941489797591 2.151929853314758
0 -0.000000000000000 24.050983364746372 5.485780545767436 8.788330738343028
0 0 1.850067196591002 0.511428337294108 0.820535933786774
6 列
-32.039999999999999
-2.058916761687573
-2.943008276653285
-27.181003032041335
-0.161268540522822
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 -0.139144811858609 1.298266818700114 0.228677309007982 0.437856328392246
-0.000000000000000 -0.000000000000000 5.250792428091452 1.397941489797591 2.151929853314758
0.000000000000002 -0.000000000000000 0 -0.917418204984064 -1.068472746647624
0.000000000000000 0.000000000000000 -0.000000000000000 0.018876836078484 0.062323742864269
6 列
-32.039999999999999
-2.058916761687573
-2.943008276653285
-13.700705710551739
0.875672672922005
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 -0.139144811858609 1.298266818700114 0.228677309007982 0.437856328392246
-0.000000000000000 -0.000000000000000 5.250792428091452 1.397941489797591 2.151929853314758
0.000000000000002 -0.000000000000000 0 -0.917418204984064 -1.068472746647624
0.000000000000000 0.000000000000000 -0.000000000000000 0 0.040338802099828
6 列
-32.039999999999999
-2.058916761687573
-2.943008276653285
-13.700705710551739
0.593766368411596
(未交换主元）Ax=b的解向量x：
-1.464947835644838
1.458145080130749
-6.004831410485430
-2.209093093570237
14.719484404672830
交换主元的高斯消元法:
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
4.930000000000000 1.210000000000000 4.480000000000000 1.100000000000000 1.000000000000000
3.530000000000000 1.460000000000000 2.920000000000000 1.210000000000000 1.000000000000000
5.050000000000000 4.040000000000000 2.510000000000000 2.010000000000000 1.000000000000000
3.540000000000000 1.040000000000000 3.470000000000000 1.020000000000000 1.000000000000000
6 列
-32.039999999999999
-20.070000000000000
-8.529999999999999
-6.300000000000000
-12.039999999999999
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 2.658015963511973 -0.749179019384265 1.117468643101482 0.424173318129989
-0.000000000000000 0.493979475484607 0.641801596351197 0.586111744583808 0.597491448118586
0 -0.139144811858609 1.298266818700114 0.228677309007982 0.437856328392246
0 0.071242873432155 1.185347776510832 0.394344355758267 0.596351197263398
6 列
-32.039999999999999
12.149486887115163
4.366374002280503
-2.058916761687573
0.892907639680731
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 2.658015963511973 -0.749179019384265 1.117468643101482 0.424173318129989
0 0 1.205428041937642 0.364392813631450 0.584982068397481
0 0 1.259047952022239 0.287175815501828 0.460061430753127
-0.000000000000000 0.000000000000000 0.781032911783379 0.378435574926643 0.518660878219538
6 列
-32.039999999999999
12.149486887115163
0.567264615543011
-1.422901745113855
2.108450332034939
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 2.658015963511973 -0.749179019384265 1.117468643101482 0.424173318129989
0 0 1.205428041937642 0.364392813631450 0.584982068397481
-0.000000000000000 0.000000000000000 0 0.142334563201356 0.139633485188985
0 0 0 -0.093425937379935 -0.150941839090899
6 列
-32.039999999999999
12.149486887115163
0.567264615543011
1.740902607064927
-2.015399453484785
1 至 5 列
8.770000000000000 2.400000000000000 5.660000000000000 1.550000000000000 1.000000000000000
0 2.658015963511973 -0.749179019384265 1.117468643101482 0.424173318129989
0 0 1.205428041937642 0.364392813631450 0.584982068397481
-0.000000000000000 0.000000000000000 0 0.142334563201356 0.139633485188985
-0.000000000000000 0.000000000000000 0 0 -0.059288842446975
6 列
-32.039999999999999
12.149486887115163
0.567264615543011
1.740902607064927
-0.872701191769358
(交换主元）Ax=b的解向量x：
-1.464947835644832
1.458145080130735
-6.004831410485433
-2.209093093570231
14.719484404672835
比较：
(未交换主元）Ax=b的解向量x：
-1.464947835644838
1.458145080130749
-6.004831410485430
-2.209093093570237
14.719484404672830
(交换主元）Ax=b的解向量x：
-1.464947835644832
1.458145080130735
-6.004831410485433
-2.209093093570231
14.719484404672835
明显我们能发现过程中的主元 存在 不同处（前者为不交换主元，后者为交换主元）的第一次消元
这说明程序的选主元模块是生效的，但结果上来看，两者间的误差也不甚明显（可能是我选的例子不好）,但也能在小数点的最后两位窥见不同。
程序优化
clc
while 1
A=input('请输入系数矩阵（方阵A）：');
b=input('请输入常数向量（列向量b）：');
[H,L,MARK]=isLegal(A,b);
if(MARK==1)
[Ax,bx]=GSelimination_NoChange(A,b,H);
[x]=GSback(Ax,bx,H);
disp("(未交换主元）Ax=b的解向量x：")
%format long
disp(x);
[Ax,bx]=GSelimination_Change(A,b,H);
[x]=GSback(Ax,bx,H);
disp("(交换主元）Ax=b的解向量x：")
%format long
disp(x);
break;
else
disp("输入不合法，请重新输入：");
end
end
function [H,L,MARK]=isLegal(A,b)
%判断输入是否合法
%合法MARK返回1，否则返回0
[HA,LA]=size(A); %返回A矩阵的行和列，分别用HA和LA表示
[Hb,Lb]=size(b); %返回b向量的行（正常来说应该为1）和列，用Hb和Lb表示
if(HA==LA)
if(HA==Hb)
if(Lb==1)
if(det(A)~=0)
MARK=1;H=HA;L=LA;
else
MARK=0;H=0;L=0;
disp("该方程组无解！");
return;
end
else
MARK=0;H=0;L=0;
disp("常数向量b输入不合法！");
return;
end
else
MARK=0;H=0;L=0;
disp("系数矩阵A和常数向量b行数不相等！");
return;
end
else
MARK=0;H=0;L=0;
disp("系数矩阵A不是方阵！");
return;
end
end
function [Ax,bx] = GSelimination_Change(A,b,n)
disp("交换主元的高斯消元法:")
%消去过程
%n表示系数矩阵(方阵)的规模
%A表示系数矩阵
%b表示常数矩阵
%消去前先交换主元
for k = 1 : n - 1
ip = k;
for i = k + 1 : n % 寻找主元的下标
if abs(A(i,k)) >= abs(A(k,k))
ip = i;
end
end
if(ip ~= k)
A([k ip],:) = A([ip k],:); %交换行
b([k ip],:) = b([ip k],:);
end
disp([A b])
if(A(k,k)~=0) %主元素不能为0,
for i=(k+1):n %循环行
c=(-1*A(i,k))/A(k,k); %倍乘因子
for j=1:n %单行循环
A(i,j)=A(i,j)+c*A(k,j); %更新矩阵元素
end
b(i)=b(i)+c*b(k); %更新常数向量
end
else
disp("主元素出现零，程序错误！"); %不更换主元，高斯消元并不总能进行
return;
end
end
disp([A b])
Ax=A;bx=b;
end
function [Ax,bx]=GSelimination_NoChange(A,b,n)
disp("不交换主元的高斯消元法：")
%消去过程
%n表示系数矩阵(方阵)的规模
%A表示系数矩阵
%b表示常数矩阵
for k=1:(n-1) %要进行n-1次消去过程
disp([A b])
if(A(k,k)~=0) %主元素不能为0,
for i=(k+1):n %循环行
c=(-1*A(i,k))/A(k,k); %倍乘因子
for j=1:n %单行循环
A(i,j)=A(i,j)+c*A(k,j); %更新矩阵元素
end
b(i)=b(i)+c*b(k); %更新常数向量
end
else
disp("主元素出现零，程序错误！"); %不更换主元，高斯消元并不总能进行
return;
end
end
disp([A b])
Ax=A;bx=b;
end
function [x]=GSback(Ax,bx,n)
%回代过程
%Ax为经过消去过程后得到的上三角矩阵
%bx为经过消去过程后得到的常数向量
%n表示系数矩阵的规模
for i=n:-1:1
if(Ax(i,i)~=0)
x(i)=bx(i);
for j=n:-1:i+1
x(i)=x(i)-Ax(i,j)*x(j);
end
x(i)=x(i)/Ax(i,i);
else
disp("主元素出现零，程序错误！");
return;
end
end
x=x';
end
328
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
