MATLAB 之线性方程组求解

2024-04-23 1109阅读

这里写目录标题

一、线性方程组求解
- 1. 线性方程组的直接解法
- - 1.1 利用左除运算符的直接解法
  - 1.2 利用矩阵的分解求解线性方程组
  - 2. 线性方程组的迭代解法
  - - 2.1 Jacobi 迭代法
    - 2.2 Gauss-Serdel 迭代法
    - 3. 求线性方程的通解
      一、线性方程组求解
      - 在 MATLAB 中，关于线性方程组的解法一般分为两类：一类是直接法，就是在没有舍入误差的情况下，通过有限步的矩阵初等运算来求得方程组的解；另一类是迭代法，就是先给定一个解的初始值，然后按照一定的迭代算法进行逐步逼近，求出更精确的近似解。
        1. 线性方程组的直接解法
        
        线性方程组的直接解法大多基于高斯消元法、主元素消元法、平方根法和追赶法等。在 MATLAB 中，这些算法已经被编成了现成的库函数或运算符，因此，只需调用相应的函数或运算符即可完成线性方程组的求解。
        1.1 利用左除运算符的直接解法
        
        线性方程组求解最简单的方法就是使用左除运算符 \，系统会自动根据输入的系数矩阵判断选用哪种方法进行求解。
        对于线性方程组 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=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >> b=[13,-9,6,0]'; >> x=A\b
        
        程序运行结果如下：
        x = -66.5556 25.6667 -18.7778 26.5556
        
        1.2 利用矩阵的分解求解线性方程组
        
        矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有 LU 分解、QR 分解、Cholesky 分解，以及 Schur 分解、Hessenberg 分解、奇异分解等。
        通过这些分解方法求解线性方程组的有点是运算速度快，可以节省存储空间。
        （1） LU 分解。矩阵的 LU 分解就是将一个矩阵表示为一个变换下三角阵和一个上三角阵的乘积形式。线性代数中已经证明，只要方阵 X X X 是非奇异的，LU 分解总是可以进行的。
        MATLAB 提供的 lu 函数用于对矩阵进行 LU 分解，其调用格式如下。
        ① [L,U]=lu(X)：产生一个上三角阵 U 和一个变换形式的下三角阵 L（行交换），使之满足 X = L U X=LU X=LU。注意，这里的矩阵 X X X 必须是方阵。
        ② [L,U,P]=lu(X)：产生一个上三角阵 U 和一个下三角阵 L 以及一个置换矩阵 P，使之满足 P X = L U PX=LU PX=LU。当然矩阵 X X X 同样必须是方阵。
        当使用第一种格式时，矩阵 L L L 往往不是一个下三角阵，但可以通过行交换成为一个下三角阵。
        例如，我们设 A = [ 1 − 1 1 5 − 4 3 2 1 1 ] A=\begin{bmatrix} 1 & -1 & 1\\ 5 & -4 & 3\\ 2 & 1 &1 \end{bmatrix} A= 152−1−41131
        则对矩阵 A A A 进行 LU 分解的程序如下：
        >> A=[1,-1,1;5,-4,3;2,1,1]; >> [L,U]=lu(A) L = 0.2000 -0.0769 1.0000 1.0000 0 0 0.4000 1.0000 0 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846
        
        为检验结果是否正确，输入如下程序：
        >> LU=L*U LU = 1 -1 1 5 -4 3 2 1 1
        
        说明结果是正确的。例中所获得的矩阵 L L L 并不是一个下三角阵，但经过各行互换之后，即可获得一个下三角阵。
        利用第二种格式对矩阵 A A A 进行 LU 分解，程序如下：
        >> [L,U,P]=lu(A) L = 1.0000 0 0 0.4000 1.0000 0 0.2000 -0.0769 1.0000 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 P = 0 1 0 0 0 1 1 0 0 >> LU=L*U %这种分解其乘积不为A LU = 5 -4 3 2 1 1 1 -1 1 >> inv(P)*L*U %考虑矩阵P后的结果 ans = 1 -1 1 5 -4 3 2 1 1
        
        实现 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.0000 0 0 0 0.5000 1.0000 0 0 0 -0.3636 0.4773 1.0000 0.5000 -1.0000 1.0000 0 U = 2.0000 1.0000 -5.0000 1.0000 0 -5.5000 2.5000 6.5000 0 0 4.0000 2.0000 0 0 0 0.4091 >> x=U\(L\b) x = -66.5556 25.6667 -18.7778 26.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 ] A=\begin{bmatrix} 1 & -1 & 1\\ 5 & -4 & 3\\ 2 & 7 &10 \end{bmatrix} A= 152−1−471310
        则对矩阵 A A A 进行 Q R QR 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.3532 0.2048 -0.3651 0.9307 -0.0228 R = -5.4772 1.2780 -6.5727 0 8.0229 8.1517 0 0 -0.5917
        
        为检验结果是否正确，输入如下程序：
        >> QR=Q*R QR = 1.0000 -1.0000 1.0000 5.0000 -4.0000 3.0000 2.0000 7.0000 10.0000
        
        说明结果是正确的。此时，我们利用第二种格式对矩阵 A A A 进行 QR 分解：
        >> [Q,R,E]=qr(A) Q = -0.0953 -0.2514 -0.9632 -0.2860 -0.9199 0.2684 -0.9535 0.3011 0.0158 R = -10.4881 -5.4347 -3.4325 0 6.0385 -4.2485 0 0 0.4105 E = 0 0 1 0 1 0 1 0 0 >> Q*R/E %验证A=Q*R*inv(E) ans = 1.0000 -1.0000 1.0000 5.0000 -4.0000 3.0000 2.0000 7.0000 10.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
        程序如下：
        >> 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)
        
        程序运行结果如下：
        x = -66.5556 25.6667 -18.7778 26.5556
        
        或采用 QR 分解的第二种格式，程序如下：
        >> [Q,R,E]=qr(A); >> x=E*(R\(Q\b)) x = -66.5556 25.6667 -18.7778 26.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 ] A=\begin{bmatrix} 2 & 1 & 1\\ 1 & 2 & -1\\ 1 & -1 &3 \end{bmatrix} A= 21112−11−13
        则对矩阵 A A A 进行 Cholesky 分解的命令如下：
        >> A=[2,1,1;1,2,-1;1,-1,3]; >> R=chol(A) R = 1.4142 0.7071 0.7071 0 1.2247 -1.2247 0 0 1.0000
        
        可以验证 R ′ R = A R'R =A R′R=A，输入如下程序：
        >> R'*R ans = 2.0000 1.0000 1.0000 1.0000 2.0000 -1.0000 1.0000 -1.0000 3.0000
        
        说明结果是正确的。此时，我们利用第二种格式对矩阵 A A A 进行 Cholesky 分解：
        >> [R,p]=chol(A) R = 1.4142 0.7071 0.7071 0 1.2247 -1.2247 0 0 1.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
        程序如下：
        >> 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) 错误使用 chol 矩阵必须为正定矩阵。
        
        命令执行时，出现错误信息，说明 A A A 为非正定矩阵。
        2. 线性方程组的迭代解法
        
        迭代法是一种不断用变量的旧值递推新值的过程，是用计算机解决问题的一种基本方法。它利用计算机运算速度快，适合做重复性操作的特点，让一组指令重复执行，在每次执行这组指令时，都从变量的原值推出它的新值。
        迭代解法非常适合求解大型稀疏矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi 迭代法、Gauss-Serdel 迭代法、超松弛迭代法和两步迭代法。
        首先，我们用一个例子说明迭代法的思想。
        为了求解线性方程组 { 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=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.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)
        
        程序运行结果如下：
        x = 0.9958 0.9579 0.7916 n = 11
        
        2.2 Gauss-Serdel 迭代法
        
        在 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
        该式即为 Gauss-Serdel 迭代公式。和 Jacobi 迭代公式相比，Gauss-Serdel 迭代用新分量代替旧分量，精度会高些。
        Gauss-Serdel 迭代法的 MATLAB 函数文件 gauseidel.m 如下：
        function [y,n]=gauseidel(A,b,x0,ep) if nargin==3 ep=1.0e-6; elseif nargin=ep x0=y; y=G*x0+f; n=n+1; end
        
        例如，我们用 Gauss-Serdel 迭代法求解下列线性方程组。设迭代初值为 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
        在程序中，我们调用函数文件 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)
        
        程序运行结果如下：
        x = 0.9958 0.9579 0.7916 n = 7
        
        由此可见，一般情况下 Gauss-Serdel 迭代比 Jacobi 迭代要收敛快一些。但这也是不绝对的，在某些情况下，Jacobi 迭代收敛而 Gauss-Serdel 迭代却可能不收敛。
        例如，我们分别用 Jacobi 迭代和 Gauss-Serdel 迭代法求解下列线性方程组，看是否收敛。 [ 1 2 − 2 1 1 1 1 2 1 ] [ x 1 x 2 x 3 ] = [ 9 7 6 ] \begin{bmatrix} 1 & 2&-2 \\ 1 & 1 & 1\\ 1&2 &1 \end{bmatrix}\begin{bmatrix}x_{1} \\x_{2} \\x_{3} \end{bmatrix}=\begin{bmatrix}9 \\7 \\6 \end{bmatrix} 111212−211 x1x2x3 = 976
        程序如下：
        >> a=[1,2,-2;1,1,1;2,2,1]; >> b=[9;7;6]; >> [x,n]=jacobi(a,b,[0;0;0]) x = -27 26 8 n = 4 >> [x,n]=gauseidel(a,b,[0;0;0]) x = NaN NaN NaN n = 1012
        
        可见对此方程，用 Jacobi 迭代收敛，而 Gauss-Serdel 迭代不收敛。因此，在使用迭代法时，要考虑算法的收敛性。
        3. 求线性方程的通解
        
        线性方程组的求解分为两类：一类是求方程组的唯一解（即特解），另一类是求方程组的无穷解（即通解）。这里对线性方程组 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') 可求得基础解系。
        （3）当方程组右端向量 b ≠ 0 b≠0 b=0 时，系数矩阵的秩 rank(4) 与其增广矩阵的秩 rank([A,b]) 是判断其是否有解的基本条件。
        ① 当 rank> b=[1;2;3]; >> [x,y]=line_solution(A,b)
        程序运行结果如下：
        方程组无解 x = [] y = []
        
        表明该方程无解。
        例如，我们求方程组的通解。 { 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 = 0 0 -8/15 3/5 y = 3/2 -3/4 3/2 7/4 1 0 0 1
        
        所以原方程的通解为 X = k 1 [ 3 / 2 3 / 2 1 0 ] + k 2 [ − 3 / 4 7 / 4 0 1 ] + [ 0 0 − 8 / 15 3 / 5 ] X=k_{1}\begin{bmatrix}3/2 \\3/2 \\1 \\0 \end{bmatrix}+k_{2}\begin{bmatrix}-3/4 \\7/4 \\0 \\1 \end{bmatrix}+\begin{bmatrix}0 \\0 \\-8/15 \\3/5 \end{bmatrix} X=k1 3/23/210 +k2 −3/47/401 + 00−8/153/5
        其中， k 1 、 k 2 k_{1}、k_{2} k1、k2 为任意常数。
        （图片来源网络，侵删）

免责声明：我们致力于保护作者版权，注重分享，被刊用文章因无法核实真实出处，未能及时与作者取得联系，或有版权异议的，请联系管理员，我们会立即处理! 部分文章是来自自研大数据AI进行生成,内容摘自(百度百科,百度知道,头条百科,中国民法典,刑法,牛津词典,新华词典,汉语词典,国家院校,科普平台)等数据,内容仅供学习参考,不准确地方联系删除处理! 图片声明：本站部分配图来自人工智能系统AI生成,觅知网授权图片,PxHere摄影无版权图库和百度，360，搜狗等多加搜索引擎自动关键词搜索配图，如有侵权的图片，请第一时间联系我们，邮箱：ciyunidc@ciyunshuju.com。本站只作为美观性配图使用,无任何非法侵犯第三方意图,一切解释权归图片著作权方,本站不承担任何责任。如有恶意碰瓷者,必当奉陪到底严惩不贷!

MATLAB 之线性方程组求解

这里写目录标题

一、线性方程组求解

1. 线性方程组的直接解法

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

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

2. 线性方程组的迭代解法

2.1 Jacobi 迭代法

2.2 Gauss-Serdel 迭代法

3. 求线性方程的通解

相关阅读

怎么把织梦的模板替换?

dedecms怎么调用特定的栏目文档?

怎么抓包一个网页?

wap怎么封装app?

目录[+]