📄 1.txt
字号:
B =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
>> [U,V,X,C,S]=gsvd(A,B)
U =
0.4082 0.7071 0.5774
-0.8165 0.0000 0.5774
0.4082 -0.7071 0.5774
V =
0.2607 -0.7950 -0.5000 0.2236
-0.4029 0.3710 -0.5000 0.6708
-0.5452 -0.0530 -0.5000 -0.6708
0.6874 0.4770 -0.5000 -0.2236
X =
0 -9.4340 -17.0587 3.4641
1.8962 8.7980 -17.0587 8.6603
3.7924 8.1620 -17.0587 13.8564
-5.6885 -7.5260 -17.0587 19.0526
C =
0 0.0000 0 0
0 0 0.0829 0
0 0 0 1.0000
S =
1.0000 0 0 0
0 1.0000 0 0
0 0 0.9966 0
0 0 0 0.0000
1.3.9 特征值问题的QZ分解
函数 qz
格式 [AA,BB,Q,Z,V] = qz(A,B) %A、B为方阵，产生上三角阵AA和BB，正交矩阵Q、Z或其列变换形式，V为特征向量阵。且满足：Q*A*Z= AA 和Q*B*Z = BB。
[AA,BB,Q,Z,V] = qz(A,B,flag) %产生由flag决定的分解结果，flag取值为：'complex'：表示复数分解（默认），取值为'real'：表示实数分解。
1.3.10 海森伯格形式的分解
如果矩阵H的第一子对角线下元素都是0，则H为海森伯格(Hessenberg)矩阵。如果矩阵是对称矩阵，则它的海森伯格形式是对角三角阵。MATLAB可以通过相似变换将矩阵变换成这种形式。
函数 hess
格式 H = hess(A) %返回矩阵A的海森伯格形式
[P,H] = hess(A) %P为酉矩阵，满足：A = PHP' 且P'P = eye(size(A))。
例1-75
>> A=[-149 -50 -154;537 180 546;-27 -9 -25];
>> [P,H]=hess(A)
P =
1.0000 0 0
0 -0.9987 0.0502
0 0.0502 0.9987
H =
-149.0000 42.2037 -156.3165
-537.6783 152.5511 -554.9272
0 0.0728 2.4489
H的第一子对角元素是H(3,1)=0。
1.4 线性方程的组的求解
我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：
若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；
若系数矩阵的秩r<n，则可能有无穷解；
线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题，通解部分属第二类问题。
1.4.1 求线性方程组的唯一解或特解（第一类问题）
这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 —— 直接法；另一类是解大型稀疏矩阵 —— 迭代法。
1．利用矩阵除法求线性方程组的特解（或一个解）
方程：AX=b
解法：X=A\b
例1-76 求方程组 的解。
解：
>>A=[5 6 0 0 0
1 5 6 0 0
0 1 5 6 0
0 0 1 5 6
0 0 0 1 5];
B=[1 0 0 0 1]';
R_A=rank(A) %求秩
X=A\B %求解
运行后结果如下
R_A =
5
X =
2.2662
-1.7218
1.0571
-0.5940
0.3188
这就是方程组的解。
或用函数rref求解：
>> C=[A,B] %由系数矩阵和常数列构成增广矩阵C
>> R=rref(C) %将C化成行最简行
R =
1.0000 0 0 0 0 2.2662
0 1.0000 0 0 0 -1.7218
0 0 1.0000 0 0 1.0571
0 0 0 1.0000 0 -0.5940
0 0 0 0 1.0000 0.3188
则R的最后一列元素就是所求之解。
例1-77 求方程组 的一个特解。
解：
>>A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
>>B=[1 4 0]';
>>X=A\B %由于系数矩阵不满秩，该解法可能存在误差。
X =[ 0 0 -0.5333 0.6000]’（一个特解近似值）。
若用rref求解，则比较精确：
>> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
B=[1 4 0]';
>> C=[A,B]; %构成增广矩阵
>> R=rref(C)
R =
1.0000 0 -1.5000 0.7500 1.2500
0 1.0000 -1.5000 -1.7500 -0.2500
0 0 0 0 0
由此得解向量X=[1.2500 – 0.2500 0 0]’（一个特解）。
2．利用矩阵的LU、QR和cholesky分解求方程组的解
（1）LU分解：
LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。
则：A*X=b 变成L*U*X=b
所以X=U\(L\b) 这样可以大大提高运算速度。
命令 [L，U]=lu (A)
例1-78 求方程组 的一个特解。
解：
>>A=[4 2 -1;3 -1 2;11 3 0];
>>B=[2 10 8]';
>>D=det(A)
>>[L,U]=lu(A)
>>X=U\(L\B)
显示结果如下：
D =
0
L =
0.3636 -0.5000 1.0000
0.2727 1.0000 0
1.0000 0 0
U =
11.0000 3.0000 0
0 -1.8182 2.0000
0 0 0.0000
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.018587e-017.
> In D:\Matlab\pujun\lx0720.m at line 4
X =
1.0e+016 *
-0.4053
1.4862
1.3511
说明 结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。
（2）Cholesky分解
若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即： 其中R为上三角阵。
方程 A*X=b 变成
所以
（3）QR分解
对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR
方程 A*X=b 变形成 QRX=b
所以 X=R\(Q\b)
上例中 [Q, R]=qr(A)
X=R\(Q\B)
说明 这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。
1.4.2 求线性齐次方程组的通解
在Matlab中，函数null用来求解零空间，即满足A?X=0的解空间，实际上是求出解空间的一组基（基础解系）。
格式 z = null % z的列向量为方程组的正交规范基，满足 。
% z的列向量是方程AX=0的有理基
例1-79 求解方程组的通解：
解：
>>A=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3];
>>format rat %指定有理式格式输出
>>B=null(A,'r') %求解空间的有理基
运行后显示结果如下：
B =
2 5/3
-2 -4/3
1 0
0 1
或通过行最简行得到基：
>> B=rref(A)
B =
1.0000 0 -2.0000 -1.6667
0 1.0000 2.0000 1.3333
0 0 0 0
即可写出其基础解系（与上面结果一致）。
写出通解：
syms k1 k2
X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解
pretty(X) %让通解表达式更加精美
运行后结果如下：
X =
[ 2*k1+5/3*k2]
[ -2*k1-4/3*k2]
[ k1]
[ k2]
% 下面是其简化形式
[2k1 + 5/3k2 ]
[ ]
[-2k1 - 4/3k2]
[ ]
[ k1 ]
[ ]
[ k2 ]
1.4.3 求非齐次线性方程组的通解
非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。
因此，步骤为：
第一步：判断AX=b是否有解，若有解则进行第二步
第二步：求AX=b的一个特解
第三步：求AX=0的通解
第四步：AX=b的通解= AX=0的通解+AX=b的一个特解。
例1-80 求解方程组
解：在Matlab中建立M文件如下：
A=[1 -2 3 -1;3 -1 5 -3;2 1 2 -2];
b=[1 2 3]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n %判断有唯一解
X=A\b
elseif R_A==R_B&R_A<n %判断有无穷解
X=A\b %求特解
C=null(A,'r') %求AX=0的基础解系
else X='equition no solve' %判断无解
end
运行后结果显示：
R_A =
2
R_B =
3
X =
equition no solve
说明 该方程组无解
例1-81 求解方程组的通解：
解法一：在Matlab编辑器中建立M文件如下：
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n
X=A\b
elseif R_A==R_B&R_A<n
X=A\b
C=null(A,'r')
else X='Equation has no solves'
end
运行后结果显示为：
R_A =
2
R_B =
2
Warning: Rank deficient, rank = 2 tol = 8.8373e-015.
> In D:\Matlab\pujun\lx0723.m at line 11
X =
0
0
-8/15
3/5
C =
3/2 -3/4
3/2 7/4
1 0
0 1
所以原方程组的通解为X=k1 +k2 +
解法二：用rref求解
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组。
运行后结果显示为：
C =
1 0 -3/2 3/4 5/4
0 1 -3/2 -7/4 -1/4
0 0 0 0 0
对应齐次方程组的基础解系为： ， 非齐次方程组的特解为： 所以，原方程组的通解为：X=k1 +k2 + 。
1.4.4 线性方程组的LQ解法
函数 symmlq
格式 x = symmlq(A,b) %求线性方程组AX=b的解X。A必须为n阶对称方阵，b为n元列向量。A可以是由afun定义并返回A*X的函数。如果收敛，将显示结果信息；如果收敛失败，将给出警告信息并显示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。
symmlq(A,b,tol) %指定误差tol，默认值是1e-6。
symmlq(A,b,tol,maxit) %maxit指定最大迭代次数
symmlq(A,b,tol,maxit,M) %M为用于对称正定矩阵的预处理因子
symmlq(A,b,tol,maxit,M1,M2) %M=M1×M2
symmlq(A,b,tol,maxit,M1,M2,x0) %x0为初始估计值，默认值为0。
[x,flag] = symmlq(A,b,…) %flag的取值为：0表示在指定迭代次数之内按要求精度收敛；1表示在指定迭代次数内不收敛；2表示M为坏条件的预处理因子；3表示两次连续迭代完全相同；4表示标量参数太小或太大；5表示预处理因子不是对称正定的。
[x,flag,relres] = symmlq(A,b,…) % relres表示相对误差norm(b-A*x)/norm(b)
[x,flag,relres,iter] = symmlq(A,b,…) %iter表示计算x的迭代次数
[x,flag,relres,iter,resvec] = symmlq(A,b,…) %resvec表示每次迭代的残差：norm(b-A*x0)
[x,flag,relres,iter,resvec,resveccg] = symmlq(A,b,…) %resveccg表示每次迭代共轭梯度残差的范数
1.4.5 双共轭梯度法解方程组
函数 bicg
格式 x = bicg(A,b) %求线性方程组AX=b的解X。A必须为n阶方阵，b为n元列向量。A可以是由afun定义并返回A*X的函数。如果收敛，将显示结果信息；如果收敛失败，将给出警告信息并显示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。
bicg(A,b,tol) %指定误差tol，默认值是1e-6。
bicg(A,b,tol,maxit) %maxit指定最大迭代次数
bicg(A,b,tol,maxit,M) %M为用于对称正定矩阵的预处理因子
bicg(A,b,tol,maxit,M1,M2) %M=M1×M2
bicg(A,b,tol,maxit,M1,M2,x0) %x0为初始估计值，默认值为0。
[x,flag] = bicg(A,b,…) %flag的取值为：0表示在指定迭代次数之内按要求精度收敛；1表示在指定迭代次数内不收敛；2表示M为坏条件的预处理因子；3表示两次连续迭代完全相同；4表示标量参数太小或太大。
[x,flag,relres] = bicg(A,b,…) % relres表示相对误差norm(b-A*x)/norm(b)
[x,flag,relres,iter] = bicg(A,b,…) %iter表示计算x的迭代次数
[x,flag,relres,iter,resvec] = bicg(A,b,…) %resvec表示每次迭代的残差：norm(b-A*x0)
例1-83 调用MATLAB6.0数据文件west0479。
>> load west0479
>> A=west0479; %将数据取为系数矩阵A。
>> b=sum (A,2); %将A的各行求和，构成一列向量。
>> X=A\b; %用“\”求AX=b的解。
>> norm(b-A*X)/norm(b) %计算解的相对误差。
ans =
1.2454e-017
>> [x,flag,relres,iter,resvec] = bicg(A,b) %用bicg函数求解。
x = （全为0，由于太长，不显示出来）
flag =
1 %表示在默认迭代次数（20次）内不收敛。
relres = %相对残差relres = norm(b-A*x)/norm(b) =norm(b)/norm(b) = 1。
1
iter = %表明解法不当，使得初始估计值0向量比后来所有迭代值都好。
0
resvec = （略） %每次迭代的残差。
>> semilogy(0:20,resvec/norm(b),'-o') %作每次迭代的相对残差图形，结果如下图。
>> xlabel('iteration number') %x轴为迭代次数。
>> ylabel('relative residual') %y轴为相对残差。
图1-1 双共轭梯度法相对误差图
1.4.6 稳定双共轭梯度方法解方程组
函数 bicgstab
格式 x =bicgstab(A,b)
bicgstab(A,b,tol)
bicgstab(A,b,tol,maxit)
bicgstab(A,b,tol,maxit,M)
bicgstab(A,b,tol,maxit,M1,M2)
bicgstab(A,b,tol,maxit,M1,M2,x0)
[x,flag] = bicgstab(A,b,…)
[x,flag,relres] = bicgstab(A,b,…)
[x,flag,relres,iter] = bicgstab(A,b,…)
[x,flag,relres,iter,resvec] = bicgstab(A,b,…)
稳定双共轭梯度法解方程组，调用方式和返回的结果形式和命令bicg一样。
例1-84
>>load west0479;
>>A=west0479;
>>b=sum(A,2);
>>[x,flag]=bicgstab(A,b)
显示结果是x的值全为0，flag=1。表示在默认误差和默认迭代次数（20次）下不收敛。
若改为：
>>[L1,U1] = luinc(A,1e-5);
>>[x1,flag1] = bicgstab(A,b,1e-6,20,L1,U1)
即指定误差，并用A的不完全LU分解因子L和U作为预处理因子M=L*U，其结果是x1的值全为0，flag=2表示预处理因子为坏条件的预处理因子。
若改为
>>[L2,U2]=luinc(A,1e-6); %稀疏矩阵的不完全LU分解。
>>[x2,flag2,relres2,iter2,resvec2]=bicgstab(A,b,1e-15,10,L2,U2)
%指定最大迭代次数为10次，预处理因子M=L*U。
>>semilogy(0:0.5:iter2,resvec2/norm(b),'-o') %每次迭代的相对残差图形，见图1-2。
>>xlabel('⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -