对于线性方程组
A
x
=
b
Ax=b
Ax=b,可以利用左除运算符 \ 求解:
x
=
A
∖
b
x=A\setminus b
x=A∖b
当系数矩阵
A
A
A 为
N
×
N
N×N
N×N 的方阵时,MATLAB 会自行用高斯消元法求解线性方程组。若右端项
b
b
b 为
N
×
1
N×1
N×1 的列向量,则
x
=
A
∖
b
x=A\setminus b
x=A∖b 可获得方程组的数值解
x
x
x(
N
×
1
N×1
N×1 的列向量)。
若右端项
b
b
b 为
N
×
M
N×M
N×M 的矩阵,则
x
=
A
∖
b
x=A\setminus b
x=A∖b 可同时获得系数矩阵
A
A
A 相同的
M
M
M 个线性方程组的数值解
x
x
x(为
N
×
M
N×M
N×M 的矩阵),即
x
(
:
,
j
)
=
A
∖
b
(
:
,
j
)
,
j
=
1
,
2
,
…
,
M
x(:,j)=A\setminus b(:,j), j=1, 2, …, M
x(:,j)=A∖b(:,j),j=1,2,…,M。
这里需要注意的是,如果矩阵
A
A
A 是奇异的或接近奇异的,则 MATLAB 会给出警告信息。
例如,我们用直接解法求解下列线性方程组。
{
2
x
1
+
x
2
−
5
x
3
+
x
4
=
13
x
1
−
5
x
2
+
7
x
4
=
−
9
2
x
2
+
x
3
−
x
4
=
6
x
1
+
6
x
2
−
x
3
−
4
x
4
=
0
\left\{\begin{matrix}2x_{1}+x_{2}-5x_{3}+x_{4}=13 \\x_{1}-5x_{2}+7x_{4}=-9 \\2x_{2}+x_{3}-x_{4}=6 \\x_{1}+6x_{2}-x_{3}-4x_{4}=0 \end{matrix}\right.
⎩⎨⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
>> A=[1,-1,1;5,-4,3;2,1,1];>>[L,U]=lu(A)
L =0.2000-0.07691.00001.0000000.40001.00000
U =5.0000-4.00003.000002.6000-0.2000000.3846
为检验结果是否正确,输入如下程序:
>> LU=L*U
LU =1-115-43211
说明结果是正确的。例中所获得的矩阵
L
L
L 并不是一个下三角阵,但经过各行互换之后,即可获得一个下三角阵。
利用第二种格式对矩阵
A
A
A 进行 LU 分解,程序如下:
>>[L,U,P]=lu(A)
L =1.0000000.40001.000000.2000-0.07691.0000
U =5.0000-4.00003.000002.6000-0.2000000.3846
P =010001100>> LU=L*U %这种分解其乘积不为A
LU =5-432111-11>>inv(P)*L*U %考虑矩阵P后的结果
ans =1-115-43211
实现 LU 分解后,线性方程组
A
x
=
b
Ax=b
Ax=b 的解
x
=
U
∖
(
L
∖
b
)
x=U\setminus (L\setminus b)
x=U∖(L∖b) 或
x
=
U
∖
(
L
∖
P
b
)
x=U\setminus (L\setminus Pb)
x=U∖(L∖Pb),这样可以大大提高运算速度。
例如,我们用 LU 分解求解下列线性方程组。
{
2
x
1
+
x
2
−
5
x
3
+
x
4
=
13
x
1
−
5
x
2
+
7
x
4
=
−
9
2
x
2
+
x
3
−
x
4
=
6
x
1
+
6
x
2
−
x
3
−
4
x
4
=
0
\left\{\begin{matrix}2x_{1}+x_{2}-5x_{3}+x_{4}=13 \\x_{1}-5x_{2}+7x_{4}=-9 \\2x_{2}+x_{3}-x_{4}=6 \\x_{1}+6x_{2}-x_{3}-4x_{4}=0 \end{matrix}\right.
⎩⎨⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
程序如下:
>> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];>> b=[13,-9,6,0]';>>[L,U]=lu(A)
L =1.00000000.50001.0000000-0.36360.47731.00000.5000-1.00001.00000
U =2.00001.0000-5.00001.00000-5.50002.50006.5000004.00002.00000000.4091>> x=U\(L\b)
x =-66.555625.6667-18.777826.5556
或采用 LU 分解的第二种格式,程序如下:
>>[L,U,P]=lu(A);>> x=U\(L\P*b);
将得到与上面相同的结果。
(2) QR 分解。对矩阵
X
X
X 进行 QR 分解,就是把
X
X
X 分解为一个正交矩阵
Q
Q
Q 和一个上三角矩阵
R
R
R 的乘积形式。QR 分解只能对方阵进行。MATLAB 的函数 qr 可用于对矩阵进行 QR 分解,其调用格式如下。
① [Q,R]=qr(X):产生一个正交矩阵
Q
Q
Q 和一个上三角阵
R
R
R,使之满足
X
=
Q
R
X=QR
X=QR。
② [Q,R,E]=qr(X):产生一个正交矩阵
Q
Q
Q、一个上三角阵
R
R
R 以及一个置换矩阵
E
E
E,使之满足
X
E
=
Q
R
XE=QR
XE=QR。
>> A=[1,-1,1;5,-4,3;2,7,10];>>[Q,R]=qr(A)
Q =-0.1826-0.0956-0.9785-0.9129-0.35320.2048-0.36510.9307-0.0228
R =-5.47721.2780-6.572708.02298.151700-0.5917
为检验结果是否正确,输入如下程序:
>> QR=Q*R
QR =1.0000-1.00001.00005.0000-4.00003.00002.00007.000010.0000
说明结果是正确的。此时,我们利用第二种格式对矩阵
A
A
A 进行 QR 分解:
>>[Q,R,E]=qr(A)
Q =-0.0953-0.2514-0.9632-0.2860-0.91990.2684-0.95350.30110.0158
R =-10.4881-5.4347-3.432506.0385-4.2485000.4105
E =001010100>> Q*R/E %验证A=Q*R*inv(E)
ans =1.0000-1.00001.00005.0000-4.00003.00002.00007.000010.0000
在实现 QR 分解后,线性方程组
A
x
=
b
Ax=b
Ax=b 的解为
x
=
R
∖
(
Q
∖
b
)
x=R\setminus (Q\setminus b)
x=R∖(Q∖b) 或
x
=
E
∖
(
R
∖
(
Q
∖
b
)
)
。
x=E\setminus (R\setminus (Q\setminus b))。
x=E∖(R∖(Q∖b))。
例如,我们用 QR 分解求解下列线性方程组。
{
2
x
1
+
x
2
−
5
x
3
+
x
4
=
13
x
1
−
5
x
2
+
7
x
4
=
−
9
2
x
2
+
x
3
−
x
4
=
6
x
1
+
6
x
2
−
x
3
−
4
x
4
=
0
\left\{\begin{matrix}2x_{1}+x_{2}-5x_{3}+x_{4}=13 \\x_{1}-5x_{2}+7x_{4}=-9 \\2x_{2}+x_{3}-x_{4}=6 \\x_{1}+6x_{2}-x_{3}-4x_{4}=0 \end{matrix}\right.
⎩⎨⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
>>[Q,R,E]=qr(A);>> x=E*(R\(Q\b))
x =-66.555625.6667-18.777826.5556
将得到与上面相同的结果。
(3) Cholesky 分解。如果矩阵
X
X
X 是对称正定的,则 Cholesky 分解将矩阵
X
X
X 分解成一个下三角阵和一个上三角阵的乘积。设上三角阵为
R
R
R,则下三角阵为其转置,即
X
=
R
′
R
X=R'R
X=R′R。MATLAB 函数 chol(X) 用于对矩阵
X
X
X 进行 Cholesky 分解,其调用格式如下。
① R=chol(X):产生一个上三角阵
R
R
R,使
R
′
R
=
X
R'R=X
R′R=X。若
X
X
X 为非对称正定,则输出一个出错信息。
② [R,p]=chol(X):这个命令格式将不输出出错信息。若
X
X
X 为对称正定的,则
p
=
0
p=0
p=0,
R
R
R 与上述格式得到的结果相同,否则
p
p
p 为一个正整数。如果
X
X
X 为满秩矩阵,则
R
R
R 为一个阶数为
q
=
p
−
1
q=p-1
q=p−1 的上三角阵,且满足
R
′
R
=
X
(
1
:
q
,
1
:
q
)
R'R =X(1:q,1:q)
R′R=X(1:q,1:q)。
>> A=[2,1,1;1,2,-1;1,-1,3];>> R=chol(A)
R =1.41420.70710.707101.2247-1.2247001.0000
可以验证
R
′
R
=
A
R'R =A
R′R=A,输入如下程序:
>> R'*R
ans =2.00001.00001.00001.00002.0000-1.00001.0000-1.00003.0000
说明结果是正确的。此时,我们利用第二种格式对矩阵
A
A
A 进行 Cholesky 分解:
>>[R,p]=chol(A)
R =1.41420.70710.707101.2247-1.2247001.0000
p =0
结果中出现
p
=
0
p=0
p=0,这表示矩阵
A
A
A 是一个正定矩阵。如果试图对一个非正定矩阵进行 Cholesky 分解,则将得出错误信息,所以,chol 函数还可以用来判定矩阵是否为正定矩阵。
实现 Cholesky 分解后,线性方程组
A
x
=
b
Ax=b
Ax=b 的解为
R
′
R
x
=
b
R'Rx=b
R′Rx=b,所以
x
=
R
∖
(
R
′
∖
b
)
。
x=R\setminus (R'\setminus b)。
x=R∖(R′∖b)。
例如,我们用 Cholesky 分解求解下列线性方程组。
{
2
x
1
+
x
2
−
5
x
3
+
x
4
=
13
x
1
−
5
x
2
+
7
x
4
=
−
9
2
x
2
+
x
3
−
x
4
=
6
x
1
+
6
x
2
−
x
3
−
4
x
4
=
0
\left\{\begin{matrix}2x_{1}+x_{2}-5x_{3}+x_{4}=13 \\x_{1}-5x_{2}+7x_{4}=-9 \\2x_{2}+x_{3}-x_{4}=6 \\x_{1}+6x_{2}-x_{3}-4x_{4}=0 \end{matrix}\right.
⎩⎨⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
为了求解线性方程组
{
10
x
1
−
x
2
=
9
−
x
1
+
10
x
2
−
2
x
3
=
7
−
2
x
2
+
10
x
3
=
6
\left\{\begin{matrix}10x_{1}-x_{2}=9 \\-x_{1}+10x_{2}-2x_{3}=7 \\-2x_{2}+10x_{3}=6 \end{matrix}\right.
⎩⎨⎧10x1−x2=9−x1+10x2−2x3=7−2x2+10x3=6
将方程改写为
{
x
1
=
10
x
2
−
2
x
3
−
7
x
2
=
10
x
1
−
9
x
3
=
1
10
(
6
+
2
x
2
)
\left\{\begin{matrix}x_{1}=10x_{2}-2x_{3}-7 \\x_{2}=10x_{1}-9 \\x_{3}=\frac{1}{10}(6+2x_{2}) \end{matrix}\right.
⎩⎨⎧x1=10x2−2x3−7x2=10x1−9x3=101(6+2x2)
这种形式的好处是将一组
x
x
x 代入右端,可以立即得到另一组
x
x
x。如果两组
x
x
x 相等,那么它就是方程组的解,不等时可以继续迭代。
例如,选取初值
x
1
=
x
2
=
x
3
=
0
x_{1}=x_{2}=x_{3}=0
x1=x2=x3=0,则经过一次迭代后, 得到
x
1
=
−
7
x_{1}=-7
x1=−7,
x
2
=
−
9
x_{2}=-9
x2=−9,
x
3
=
0.6
x_{3}=0.6
x3=0.6,然后再继续迭代。可以构造方程的迭代公式为
{
x
1
(
k
+
1
)
=
10
x
2
(
k
)
−
2
x
3
(
k
)
−
7
x
2
(
k
+
1
)
=
10
x
1
(
k
)
−
9
x
3
(
k
+
1
)
=
1
10
(
6
+
2
x
2
(
k
)
)
\left\{\begin{matrix}x_{1}^{(k+1)}=10x_{2}^{(k)}-2x_{3}^{(k)}-7 \\x_{2}^{(k+1)}=10x_{1}^{(k)}-9 \\x_{3}^{(k+1)}=\frac{1}{10}(6+2x_{2}^{(k)}) \end{matrix}\right.
⎩⎨⎧x1(k+1)=10x2(k)−2x3(k)−7x2(k+1)=10x1(k)−9x3(k+1)=101(6+2x2(k))
2.1 Jacobi 迭代法
对于线性方程组
A
x
=
b
Ax=b
Ax=b,如果
A
A
A 是非奇异方阵,且
a
i
i
≠
0
(
i
=
1
,
2
,
…
,
n
)
a_{ii}\ne 0(i=1,2,…,n)
aii=0(i=1,2,…,n),则可将
A
A
A 分解为
A
=
D
−
L
−
U
A=D-L-U
A=D−L−U,其中
D
D
D 为对角阵,其元素为
A
A
A 的对角元素,
L
L
L 与
U
U
U 为
A
A
A 的下三角阵取反和上三角阵取反:
L
=
−
[
0
a
2
,
1
0
⋮
⋱
⋱
a
n
,
1
⋯
a
n
,
n
−
1
0
]
L=-\begin{bmatrix} 0 & & & \\ a_{2,1} & 0 & & \\ \vdots & \ddots & \ddots & \\ a_{n,1} & \cdots & a_{n,n-1} &0 \end{bmatrix}
L=−0a2,1⋮an,10⋱⋯⋱an,n−10
U
=
−
[
0
a
1
,
2
⋯
a
1
,
n
0
⋱
⋮
⋱
a
n
−
1
,
n
0
]
U=-\begin{bmatrix} 0 & a_{1,2} & \cdots &a_{1,n} \\ & 0 & \ddots& \vdots\\ & & \ddots &a_{n-1,n} \\ & & &0 \end{bmatrix}
U=−0a1,20⋯⋱⋱a1,n⋮an−1,n0
于是
A
x
=
b
Ax=b
Ax=b 转化为
x
=
D
−
1
(
L
+
U
)
x
+
D
−
1
b
x=D^{-1}(L+U)x+D^{-1}b
x=D−1(L+U)x+D−1b
与之对应的迭代公式为
x
k
+
1
=
D
−
1
(
L
+
U
)
x
k
+
D
−
1
b
x^{k+1}=D^{-1}(L+U)x^{k}+D^{-1}b
xk+1=D−1(L+U)xk+D−1b
这就是 Jacobi 迭代公式。如果序列
{
x
k
+
1
}
\left \{x^{k+1}\right \}
{xk+1} 收敛于
x
x
x,则
x
x
x 必是方程
A
x
=
b
Ax=b
Ax=b 的解。
Jacobi 迭代法的 MATLAB 函数文件 jacobi.m 如下:
function [y,n]=jacobi(A,b,x0,ep)if nargin==3
ep=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A));%求A的对角矩阵
L=-tril(A,-1);%求A的下三角阵
U=-triu(A,1);%求A的上三角阵
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;%迭代次数
whilenorm(y-x0)>=ep
x0=y;
y=B*x0+f;
n=n+1;
end
例如,我们用 Jacobi 迭代法求解下列线性方程组。设迭代初值为 0,迭代精度为
1
0
−
6
10^{-6}
10−6。
{
10
x
1
−
x
2
=
9
−
x
1
+
10
x
2
−
2
x
3
=
7
−
2
x
2
+
10
x
3
=
6
\left\{\begin{matrix}10x_{1}-x_{2}=9 \\-x_{1}+10x_{2}-2x_{3}=7 \\-2x_{2}+10x_{3}=6 \end{matrix}\right.
⎩⎨⎧10x1−x2=9−x1+10x2−2x3=7−2x2+10x3=6
在 Jacobi 迭代过程中,计算
x
i
(
k
+
1
)
x^{(k+1)}_{i}
xi(k+1) 时,
x
1
(
k
+
1
)
,
…
,
x
i
−
1
(
k
+
1
)
x^{(k+1)}_{1},…,x^{(k+1)}_{i-1}
x1(k+1),…,xi−1(k+1) 已经得到,不必再用
x
1
(
k
)
,
…
,
x
i
−
1
(
k
)
x^{(k)}_{1},…,x^{(k)}_{i-1}
x1(k),…,xi−1(k),即原来的迭代公式
D
x
(
k
+
1
)
=
(
L
+
U
)
X
(
k
)
+
b
Dx^{(k+1)}=(L+U)X^{(k)}+b
Dx(k+1)=(L+U)X(k)+b 可以改进为
D
x
(
k
+
1
)
=
L
(
k
+
1
)
+
X
(
k
)
+
b
Dx^{(k+1)}=L^{(k+1)}+X^{(k)}+b
Dx(k+1)=L(k+1)+X(k)+b,于是得到
x
(
k
+
1
)
=
(
D
−
L
)
−
1
U
x
(
k
)
+
(
D
−
L
)
−
1
b
x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b
x(k+1)=(D−L)−1Ux(k)+(D−L)−1b
线性方程组的求解分为两类:一类是求方程组的唯一解(即特解),另一类是求方程组的无穷解(即通解)。这里对线性方程组
A
x
=
b
Ax=b
Ax=b 的求解理论进行归纳。
(1) 当系数矩阵
A
A
A 是一个满秩方阵时,方程
A
x
=
b
Ax=b
Ax=b 称为恰定方程,方程有唯一 解
x
=
A
−
1
b
x=A^{-1}b
x=A−1b,这是最基本的一种情况。一般用
x
=
A
∖
b
x=A\setminus b
x=A∖b 求解速度更快。
(2) 当方程组右端向量
b
=
0
b=0
b=0 时,方程称为齐次方程组。齐次方程组总有零解,因此称解
x
=
0
x=0
x=0 为平凡解。当系数矩阵
A
A
A 的秩小于
n
n
n(
n
n
n 为方程组中未知变量的个数)时,齐次方程组有无穷多个非平凡解,其通解中包含
n
n
n-rank(4) 个线性无关的解向量,用 MATLAB 的函数 null(A,'r') 可求得基础解系。
① 当 rank<(A)=rank([4,b])=n 时,方程组有唯一解:
x
r
=
A
∖
b
xr=A\setminus b
xr=A∖b 或
x
=
p
i
n
v
(
A
)
∗
b
x=pinv(A)*b
x=pinv(A)∗b。
② 当 rank(4)=rank([A,b])<n 时,方程组有无穷多个解,其通解 = 方程组的一个特解 + 对应的齐次方程组
A
x
=
0
Ax=0
Ax=0 的通解。可以用
A
∖
b
A\setminus b
A∖b 求得方程组的一个特解,用 null(A,'r') 求得该方程组所对应的齐次方程组的基础解系,基础解系中包含
n
n
n-rank(4) 个线性无关的解向量。
③ 当 rank(A)<rank([A,b]) 时,方程组无解。
下面,我们写一个求解线性方程组的函数文件 line_solution.m。
function [x,y]=line_solution(A,b)[m,n]=size(A);
y=[];ifnorm(b)>0%非齐次方程组
ifrank(A)==rank([A,b])ifrank(A)==n %有唯一解
disp('原方程组有唯一解x');
x=A\b;else%方程组有无穷多个解,基础解系
disp('原方程组有无穷个解,特解为x,其齐次方程组的基础解系为y');
x=A\b;
y=null(A,'r');
end
elsedisp('方程组无解');%方程组无解
x=[];
end
else%齐次方程组
disp('原方程组有零解x');
x=zeros(n,1);%0解
ifrank(A)<n
disp('方程组有无穷个解,基础解系为y');%非0解
y=null(A,'r');
end
end
例如,我们求解方程组。
{
x
1
−
2
x
2
+
3
x
3
−
x
4
=
1
3
x
1
−
x
2
+
5
x
3
−
3
x
4
=
2
2
x
1
+
x
2
+
2
x
3
−
2
x
4
=
3
\left\{\begin{matrix}x_{1}-2x_{2}+3x_{3}-x_{4}=1 \\3x_{1}-x_{2}+5x_{3}-3x_{4}=2 \\2x_{1}+x_{2}+2x_{3}-2x_{4}=3 \end{matrix}\right.
⎩⎨⎧x1−2x2+3x3−x4=13x1−x2+5x3−3x4=22x1+x2+2x3−2x4=3
例如,我们求方程组的通解。
{
x
1
+
x
2
−
3
x
3
−
x
4
=
1
3
x
1
−
x
2
−
3
x
3
+
4
x
4
=
4
x
1
+
5
x
2
−
9
x
3
−
8
x
4
=
0
\left\{\begin{matrix}x_{1}+x_{2}-3x_{3}-x_{4}=1 \\3x_{1}-x_{2}-3x_{3}+4x_{4}=4 \\x_{1}+5x_{2}-9x_{3}-8x_{4}=0 \end{matrix}\right.
⎩⎨⎧x1+x2−3x3−x4=13x1−x2−3x3+4x4=4x1+5x2−9x3−8x4=0
程序如下:
>> format rat %指定有理式格式输出
>> A=[1,1,-3,-1;3,-1,-3,4;1,5,-9,-8];>> b=[1;4;0];>>[x,y]=line_solution(A,b)>> x,y
>> format short%恢复默认的短格式输出
程序运行结果如下:
警告: 秩亏,秩 =2,tol =8.837264e-15。
> 位置:line_solution(第 11 行)
x =00-8/153/5
y =3/2-3/43/27/41001