2.1 例题解答
例 2.1 用Gauss消元法解方程组:
解:
直接建立求解该方程组的M文件Gauss.m如下:
%求解例题2.1
%高斯法求解线性方程组Ax=b
%A为输入矩阵系数,b为方程组右端系数
%方程组的解保存在x变量中
%先输入方程系数
A=[1 2 3;2 7 5;1 49];
b=[1 6 -3]';
[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 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);
直接运行上面的M文件或在Matlab命令窗口中直接输入Gauss即可得出结果.
在Matlab命令窗口中输入Gauss得出结果如下:
>> Gauss
x=
2.0000
1.0000
-1.0000
扩展:
Matlab求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成: ,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一组解):
运用求逆思想 : 或 ;
左除法,原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过 分解法进行的(即先将矩阵A作 分解,再回代计算): ;
用符号矩阵法计算,这种计算方法最接近精确值,但计算速度最慢:
通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以这样实现之: ; .
上面四种常用的办法示例如下:
>> A=[1 23;2 7 5;1 4 9] %上面示例方程组系数
A =
1 2 3
2 7 5
1 4 9
>> b=[1 6-3]' %方程组右端的系数
b =
1
6
-3
>>x1_1=inv(A)*b,x1_2=A^(-1)*b %方法一,求逆思想
x1_1 =
2.0000
1.0000
-1.0000
x1_2 =
2.0000
1.0000
-1.0000
>> x2=A\b %方法二,左除思想
x2 =
2
1
-1
>>x3=sym(A)\sym(b) %方法三,符号法
x3 =
2
1
-1
>>C=[A,b],rref(C) %方法四, 行简化阶梯形思想,最后输出结果的一列为解
C =
1 2 3 1
2 7 5 6
1 4 9 -3
ans =
1 0 0 2
0 1 0 1
0 0 1 -1
例 2.2 用选列主元的Gauss消元法解如下方程:
解:
直接建立求解该方程的M文件Gauss_line.m,求解程序编制如下:
%求解例题2.2
%高斯列主元消元法求解线性方程组Ax=b
%A为输入矩阵系数,b为方程组右端系数
%方程组的解保存在x变量中
format long;%设置为长格式显示,显示15位小数
A=[0.00001,1;2,1]
b=[1.00001,3]'
[m,n]=size(A);
%先检查系数正确性
if m~=n
error('矩阵A的行数和列数必须相同');
return;
end
if m~=size(b)
error('b的大小必须和A的行数或A的列数相同');
return;
end
%再检查方程是否存在唯一解
ifrank(A)~=rank([A,b])
error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');
return;
end
c=n+1;
A(:,c)=b; %(增广)
for k=1:n-1
[r,m]=max(abs(A(k:n,k))); %选主元
m=m+k-1; %修正操作行的值
if(A(m,k)~=0)
if(m~=k)
A([k m],:)=A([m k],:); %换行
end
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
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);
format short;%设置为默认格式显示,显示5位
运行得到输出结果如下:
A =
0.000010000000000 1.000000000000000
2.000000000000000 1.000000000000000
b =
1.000010000000000
3.000000000000000
X=
1.000000000000000
1.000000000000000
例 2.3 用 消元法思想求
的逆阵.
解:
解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下:
%Gauss-Jordan法求例2.3
clc;
A=[1 -1 0;2 2 3;-12 1];
A1=A;%先保存原来的方阵A
[n,m]=size(A);
if n~=m
error('A必须为方阵');
return;
end
A(:,n+1:2*n)=eye(n);%构造增广矩阵
for k=1:n
[l,m]=max(abs(A(k:n,k)));%按列选主元
if A(k+m-1,k)==0
error('找到列最大的元素为零,错误');
return;
end
if m~=1 %交换
Temp=A(k,:);
A(k,:)=A(k+m-1,:);
A(k+m-1,:)=Temp;
end
for i=1:n
if i~=k
A(i,:)=A(i,:)-A(k,:)*A(i,k)/A(k,k);
end
end
end
for i=n:(-1):1
A(i,:)=A(i,:)/A(i,i);
end
A(:,1:n)=[];
disp('A=');
disp(A1);
disp('用Gauss-Jandan算得矩阵A的逆矩阵为:');
disp('inv(A)=');
disp(A);
clear Temp i k l mn;%清除临时变量
在Matlab命令窗口中输入Gauss_Jordan回车后得到结果如下:
A=
1 -1 0
2 2 3
-1 2 1
用Gauss-Jandan算得矩阵A的逆矩阵为:
inv(A)=
-4 1 -3
-5 1 -3
6 -1 4
解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助 命令求解,非常方便:
输入矩阵:
>> A=[1 -10;2 2 3;-1 2 1]
A =
1 -1 0
2 2 3
-1 2 1
>>C=[A,eye(length(A))]
C =
1 -1 0 1 0 0
2 2 3 0 1 0
-1 2 1 0 0 1
>>invA=rref(C)
invA =
1 0 0 -4 1 -3
0 1 0 -5 1 -3
0 0 1 6 -1 4
后三列即为其逆阵,检验其正确性:
>>c=invA(:,4:6)
c =
-4 1 -3
-5 1 -3
6 -1 4
>> d=c*A
d =
1 0 0
0 1 0
0 0 1
可见求解正确.
例2.4用分解LU的方法求解方程组
解:
解线性方程组中 分解的L,U可以实现矩阵A的三角分解,使得A=L*U. L,U应该是下三角和上三角矩阵的,这样才利于回代求根.但是MATLAB中的LU分解与解线性方程组中的L,U不一样.
MATLAB的LU分解命令调用格式为:
[L,U]=lu(A)
MATLAB计算出来的L是"准下三角"(交换L的行后才能成为真正的下三角阵),U为上三角矩阵,但它们还是满足A=L*U的.
先录入矩阵系数.
>> A=[2 4 26;4 9 6 15;2 6 9 18;6 15 18 40]
A =
2 4 2 6
4 9 6 15
2 6 9 18
6 15 18 40
>> b=[9 2322 47]'
b =
9
23
22
47
将A作 分解,方法是使用矩阵分解的LU命令即可:
>>[L,D]=lu(A)
L =
0.3333 1.0000 -0.6667 1.0000
0.6667 1.0000 0 0
0.3333 -1.0000 1.0000 0
1.0000 0 0 0
U =
6.0000 15.0000 18.0000 40.0000
0 -1.0000 -6.0000 -11.6667
0 0 -3.0000 -7.0000
0 0 0 -0.3333
再检验其正确性:
>>C=L*U
C =
2 4 2 6
4 9 6 15
2 6 9 18
6 15 18 40
解方程组
>>y=L\b
y =
47.0000
-8.3333
-2.0000
0.3333
解方程组 得到方程组的最终解:
>>x=U\y
x =
0.5000
2.0000
3.0000
-1.0000
故方程组的最终解为:
例 2.5 解方程组
试用平方根法,改进的平方根法和追赶法分别解之.
解:
先输入矩阵:
>>A=[6 1 0;1 4 1;0 1 14]
A =
6 1 0
1 4 1
0 1 14
>>b=[6 24 322]'
b =
6
24
322
平方根法:
先对A矩阵作 分解:
>>L=chol(A)
L =
2.4495 0.4082 0
0 1.9579 0.5108
0 0 3.7066
检验其正确性
>>L'*L
ans =
6.0000 1.0000 0
1.0000 4.0000 1.0000
0 1.0000 14.0000
将L转化为下三角矩阵:
>>L=L'
L =
2.4495 0 0
0.4082 1.9579 0
0 0.5108 3.7066
解方程组 :
>>y=L\b
y =
2.4495
11.7473
85.2526
再解方程组 ,得到最终解:
>>x=L'\y
x =
1.0000
-0.0000
23.0000
故平方根法解上述方程的结果为:
改进的平方根法:
先对矩阵A作LDL分解:
>> [L,D]=ldl(A)
L =
1.0000 0 0
0.1667 1.0000 0
0 0.2609 1.0000
D =
6.0000 0 0
0 3.8333 0
0 0 13.7391
检验其分解正确性:
>>L*D*L'
ans =
6 1 0
1 4 1
0 1 14
解方程组 :
>> y=L\b
y =
6
23
316
解方程组 :
>> x=(D*L')\y
x =
1
0
23
故改进的平方根法解上述方程的结果亦为:
追赶法:
编制追赶法求解该方程的程序如下:
%pursue.m
%三对角线性方程组的追赶法解方程组例2.5
%输入矩阵
clc;
A=[6 1 0;1 4 1;0 114]
f=[6 24 322]
[n,m]=size(A);
%分别取对角元素
a=zeros(1,n);
a(2:n)=diag(A,-1);
c=diag(A,1);
%此处用变量d存储A主对角线上的元素,因已用变量b存储方程右边的系数
b=diag(A);
if b(1)==0
error('主对角元素不能为0');
return;
end
%初始计算,式(2.31)
alpha(1)=b(1);
beta(1)=c(1)/b(1);
%按照公式(2.31)计算
for i=2:n-1
alpha(i)=b(i)-a(i)*beta(i-1);
if alpha(i)==0
error('错误:在解方程过程中α为0');
return;
end
beta(i)=c(i)/alpha(i);
end
%对最后一行作计算
alpha(n)=b(n)-a(n)*beta(n-1);
if alpha(n)==0
error('错误:在解方程过程中最后一个α为0');
return;
end
%以下按照公式(2.32)计算,解Ly=f
y(1)=f(1)/b(1);
for i=2:n
y(i)=(f(i)-a(i)*y(i-1))/alpha(i);
end
%以下按照公式(2.33)计算,解Ux=y
X(n)=y(n);
for i=n-1:-1:1
X(i)=y(i)-beta(i)*X(i+1);
end
disp('X=');
disp(X);
在Matlab命令窗口输入pursue计算结果如下:
>>pursue
A =
6 1 0
1 4 1
0 1 14
f =
6 24 322
X=
1 0 23
其中A为系数矩阵,f为矩阵右端的系数,最后计算结果为X.
由以上计算可知追赶法解该方程的结果亦为:
例 2.6 ,求 .
解:
输入矩阵:
>>x=[1 0.5 -0.3]
x =
1.0000 0.5000 -0.3000
计算其1-范数
>>norm_1=norm(x,1)
norm_1 =
1.8000
计算其2-范数
>>norm_2=norm(x)
norm_2 =
1.1576
计算其无穷大范数
>>norm_inf=norm(x,inf)
norm_inf =
1
还可以计算其无穷小范数(即各分量绝对值中的最小值)
>>norm_minusInf=norm(x,-inf)
norm_minusInf =
0.3000
例 2.7
求 .
解:
先输入矩阵:
>>A=[-2 1 0;1 -2 1;0 1 -2]
A =
-2 1 0
1 -2 1
0 1 -2
求A的1-范数(列和范数):
>>norm_1=norm(A,1)
norm_1 =
4
求解A的无穷大范数(行和范数):
>>norm_inf=norm(A,inf)
norm_inf =
4
求A的2-范数( 最大特征值):
>>norm_2=norm(A,2)
norm_2 =
3.4142
还可以求解A的F-范数:
>>norm_F=norm(A,'fro')
norm_F =
4
谱半径可以通过按其概念进行计算:对其特征值的绝对值取最大值即可:
>>R_A=max(abs(eig(A)))
R_A =
3.4142
例 2.8 n阶Hilbert矩阵
考查其条件数.
解:
上述矩阵为抽象矩阵,而计算机只能进行有限次迭代,我们从 考查其条件数,为此编制如下程序Hilb_cond10.m:
% Hilb_cond10.m
%考查从3阶至10的 矩阵2—条件数
for i=3:10
h=hilb(i);
condA(i)=cond(h,2);
end
disp('n cond');
for i=3:10
s=sprintf('%d %f',i,condA(i));
disp(s);
end
运行后得到如下结果:
n cond
3 524.056778
4 15513.738739
5 476607.250242
6 14951058.641005
7 475367356.370393
8 15257575270.772364
9 493153986466.270940
10 16025391750078.617000
结果中左边为 的阶数,右边为对应的条件数 ,从以上结果也可分析可知:当n的阶数稍高时, 矩阵即出现严重的病态
.
例 2.9 求
的条件数
解:
建立2阶 矩阵:
>>H=hilb(2)
H =
1.0000 0.5000
0.5000 0.3333
求其2-条件数
>>cond_2=cond(H,2)
cond_2 =
19.2815
求其1-条件数
>>cond_1=cond(H,1)
cond_1 =
27.0000
求其无穷大条件数
>>cond_inf=cond(H,inf)
cond_inf =
27.0000
还可求其F-条件数
>>cond_f=cond(H,'fro')
cond_f =
19.3333
