Matlab编程资源库（12）线性方程组求解

一、直接解法

1 ．利用左除运算符的直接解法

     对于线性方程组 Ax = b ，可以利用 左除运算符 “ \ ” 求解：

x=A\b

2 ．利用矩阵的分解求解线性方程组

     矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU 分解、 QR 分解、Cholesky分解，以及 Schur 分解、 Hessenberg 分解、奇异分解等。

 (1) LU 分解

      矩阵的 LU 分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A 是非奇异的， LU 分解总是可以进行的。MATLAB提供的 lu 函数用于对矩阵进行 LU 分解，其调用格式为：

[L,U]=lu(X) ：产生一个上三角阵 U 和一个变换形式的下三角阵L( 行交换 ) ，使之满足 X=LU 。

 注意， 这里的矩阵X必须是方阵 。

[L,U,P]=lu(X) ：产生一个上三角阵 U 和一个下三角阵 L 以及一个置换矩阵P ，使满足 PX=LU 。  当然矩阵 X 同样必须是方阵。

       实现LU 分解后，线性方程组 Ax=b 的解 x=U\(L\b) 或x=U\(L\Pb)，这样可以大大提高运算速度。

  (2) QR 分解

       对矩阵 X 进行 QR 分解，就是把 X 分解为一个 正交矩阵Q 和一个上 三角矩阵R 的乘积形式。 QR 分解只能对 方阵 进行。

       MATLAB的函数 qr 可用于对矩阵进行 QR 分解，其调用格式为：

[Q,R]=qr(X) ：产生一个一个正交矩阵 Q 和一个上三角矩阵R，使之满足 X=QR 。

[Q,R,E]=qr(X) ：产生一个一个正交矩阵 Q 、一个上三角矩阵R以及一个 置换矩阵E ，使之满足 XE=QR 。

      实现 QR 分解后，线性方程组 Ax=b 的解 x=R\(Q\b) 或x=E(R\(Q\b))。

 (3) Cholesky 分解

      如果矩阵 X 是 对称正定 的，则 Cholesky 分解将矩阵 X 分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R ，则下三角矩阵为其转置，即X=R'R 。

      MATLAB 函数 chol(X)用于对矩阵X 进行 Cholesky 分解，其调用格式为：

R=chol(X) ：产生一个上三角阵 R ，使 R ' R = X 。

      若 X 为非对称正定，则输出一个出错信息。

      [R,p]=chol(X) ：这个命令格式将不输出出错信息。当 X 为对称正定的，则p=0 ， R 与上述格式得到的结果相同；否则 p为一个正整数。如果X 为满秩矩阵，则 R 为一个阶数为q=p-1的上三角阵，且满足 R'R=X(1:q,1:q) 。

      实现Cholesky 分解后，线性方程组 Ax=b 变成 R‘Rx=b ，所以x=R\(R’\b)。

例 7-1 用直接解法求解下列线性方程组。

命令如下：

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
x=A\b

运行结果 ：

例 7-2 用 LU 分解求解例 7-1 中的线性方程组。

命令如下：

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);
x=U\(L\b)

运行结果 ：

或采用 LU 分解的第 2 种格式，命令如下：

[L,U ,P]=lu(A);
x=U\(L\P*b)

运行结果 ：

例 7-3 用 QR 分解求解例 7-1 中的线性方程组。

命令如下：

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[Q,R]=qr(A);
x=R\(Q\b)

运行结果 ：

或采用 QR 分解的第 2 种格式，命令如下：

[Q,R,E]=qr(A);
x=E*(R\(Q\b))

运行结果 ：

例 7-4 用 Cholesky 分解求解例 7-1 中的线性方程组。

命令如下：

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
R=chol(A)
%{
??? Error using ==> chol
Matrix must be positive definite
命令执行时，出现错误信息，说明A为非正定矩阵。
%}

运行结果 ：

二、迭代解法

      迭代解法非常适合求解 大型系数矩阵 的方程组。在数值分析 中，迭代解法主要包括 Jacobi 迭代法、 Gauss-Serdel 迭代 法、超松弛迭代法和两步迭代法。

1 ． Jacobi 迭代法

      对于线性方程组 Ax=b ，如果 A 为非奇异方阵，即 aii ≠ 0(i=1,2,…,n) ，则可将 A 分解为 A=D-L-U ，其中 D 为 对角阵 ，其元素为 A 的对角元素， L 与 U 为 A 的下三角阵和上 三角阵，于是 Ax=b 化为：

x=D -1 (L+U)x+D -1 b

       与之对应的迭代公式为：

x(k+1)=D -1 (L+U)x(k)+D -1 b

       这就是 Jacobi 迭代公式, 如果序列 {x(k+1)} 收敛于 x ，则 x 必是 方程 Ax=b 的解。

2 ． Gauss-Serdel 迭代法

       在 Jacobi 迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx (k+1) =(L+U)x (k) +b 可以改进为Dx(k+1) =Lx (k+1 )+Ux (k) +b ，于是得到：

x (k+1) =(D-L) -1 Ux (k) +(D-L) -1 b

       该式即为 Gauss-Serdel 迭代公式。和 Jacobi 迭代相比， GaussSerdel迭代用新分量代替旧分量，精度会高些。

Jacobi 迭代法的 MATLAB 函数文件 Jacobi.m 如下：

function [y,n]=jacobi(A,b,x0,eps)
if nargin==3
eps=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; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=B*x0+f;
n=n+1;
end

例 7-5 用 Jacobi 迭代法求解下列线性方程组。设迭代初值为 0 ，迭代精度为 10 -6 。

在命令中调用函数文件 Jacobi.m ，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];
b=[9,7,6]';
[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)

运行结果 ：

Gauss-Serdel 迭代法的 MATLAB 函数文件 gauseidel.m 如下：

function [y,n]=gauseidel(A,b,x0,eps)
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end 
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=G*x0+f;
n=n+1;
end

例 7-6 用 Gauss-Serdel 迭代法求解下列线性方程组。设迭代 初值为 0 ，迭代精度为 10 -6 。

在命令中调用函数文件 gauseidel.m ，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];
b=[9,7,6]';
[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)

运行结果 ：

例 7-7 分别用 Jacobi 迭代和 Gauss-Serdel 迭代法求解下列线性 方程组，看是否收敛。

命令如下：

a=[1,2,-2;1,1,1;2,2,1];
b=[9;7;6];
[x,n]=jacobi(a,b,[0;0;0])
[x,n]=gauseidel(a,b,[0;0;0])

运行结果 ：

 结语   

不愿迈出前行的脚步，就无法到达最美的远方

不肯跳出眼前的安逸，就无法感受生活的多彩

！！！