沧州市高速公路建设管理局网站,下陆区建设局网站,wordpress的主题切换不成功,宝安区哪一个街道最富裕这里写目录标题 一、线性方程组求解1. 线性方程组的直接解法1.1 利用左除运算符的直接解法1.2 利用矩阵的分解求解线性方程组 2. 线性方程组的迭代解法2.1 Jacobi 迭代法2.2 Gauss-Serdel 迭代法 3. 求线性方程的通解 一、线性方程组求解
在 MATLAB 中#xff0c;关于线性方程… 这里写目录标题 一、线性方程组求解1. 线性方程组的直接解法1.1 利用左除运算符的直接解法1.2 利用矩阵的分解求解线性方程组 2. 线性方程组的迭代解法2.1 Jacobi 迭代法2.2 Gauss-Serdel 迭代法 3. 求线性方程的通解 一、线性方程组求解
在 MATLAB 中关于线性方程组的解法一般分为两类一类是直接法就是在没有舍入误差的情况下通过有限步的矩阵初等运算来求得方程组的解另一类是迭代法就是先给定一个解的初始值然后按照一定的迭代算法进行逐步逼近求出更精确的近似解。
1. 线性方程组的直接解法
线性方程组的直接解法大多基于高斯消元法、主元素消元法、平方根法和追赶法等。在 MATLAB 中这些算法已经被编成了现成的库函数或运算符因此只需调用相应的函数或运算符即可完成线性方程组的求解。
1.1 利用左除运算符的直接解法
线性方程组求解最简单的方法就是使用左除运算符 \系统会自动根据输入的系数矩阵判断选用哪种方法进行求解。对于线性方程组 A x b Axb Axb可以利用左除运算符 \ 求解 x A ∖ b xA\setminus b xA∖b当系数矩阵 A A A 为 N × N N×N N×N 的方阵时MATLAB 会自行用高斯消元法求解线性方程组。若右端项 b b b 为 N × 1 N×1 N×1 的列向量则 x A ∖ b xA\setminus b xA∖b 可获得方程组的数值解 x x x N × 1 N×1 N×1 的列向量。若右端项 b b b 为 N × M N×M N×M 的矩阵则 x A ∖ b xA\setminus b xA∖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), j1, 2, …, M x(:,j)A∖b(:,j),j1,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. ⎩ ⎨ ⎧2x1x2−5x3x413x1−5x27x4−92x2x3−x46x16x2−x3−4x40程序如下 A[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b[13,-9,6,0];xA\b程序运行结果如下
x -66.555625.6667-18.777826.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 XLU XLU。注意这里的矩阵 X X X 必须是方阵。② [L,U,P]lu(X)产生一个上三角阵 U 和一个下三角阵 L 以及一个置换矩阵 P使之满足 P X L U PXLU PXLU。当然矩阵 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.00001.0000 0 00.4000 1.0000 0U 5.0000 -4.0000 3.00000 2.6000 -0.20000 0 0.3846
为检验结果是否正确输入如下程序 LUL*ULU 1 -1 15 -4 32 1 1
说明结果是正确的。例中所获得的矩阵 L L L 并不是一个下三角阵但经过各行互换之后即可获得一个下三角阵。利用第二种格式对矩阵 A A A 进行 LU 分解程序如下 [L,U,P]lu(A)L 1.0000 0 00.4000 1.0000 00.2000 -0.0769 1.0000U 5.0000 -4.0000 3.00000 2.6000 -0.20000 0 0.3846P 0 1 00 0 11 0 0 LUL*U %这种分解其乘积不为ALU 5 -4 32 1 11 -1 1 inv(P)*L*U %考虑矩阵P后的结果ans 1 -1 15 -4 32 1 1
实现 LU 分解后线性方程组 A x b Axb Axb 的解 x U ∖ ( L ∖ b ) xU\setminus (L\setminus b) xU∖(L∖b) 或 x U ∖ ( L ∖ P b ) xU\setminus (L\setminus Pb) xU∖(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. ⎩ ⎨ ⎧2x1x2−5x3x413x1−5x27x4−92x2x3−x46x16x2−x3−4x40程序如下 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 00.5000 1.0000 0 00 -0.3636 0.4773 1.00000.5000 -1.0000 1.0000 0U 2.0000 1.0000 -5.0000 1.00000 -5.5000 2.5000 6.50000 0 4.0000 2.00000 0 0 0.4091 xU\(L\b)x -66.555625.6667-18.777826.5556
或采用 LU 分解的第二种格式程序如下 [L,U,P]lu(A);xU\(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 XQR XQR。② [Q,R,E]qr(X)产生一个正交矩阵 Q Q Q、一个上三角阵 R R R 以及一个置换矩阵 E E E使之满足 X E Q R XEQR XEQR。例如我们设 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.0228R -5.4772 1.2780 -6.57270 8.0229 8.15170 0 -0.5917
为检验结果是否正确输入如下程序 QRQ*RQR 1.0000 -1.0000 1.00005.0000 -4.0000 3.00002.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.0158R -10.4881 -5.4347 -3.43250 6.0385 -4.24850 0 0.4105E 0 0 10 1 01 0 0 Q*R/E %验证AQ*R*inv(E)ans 1.0000 -1.0000 1.00005.0000 -4.0000 3.00002.0000 7.0000 10.0000
在实现 QR 分解后线性方程组 A x b Axb Axb 的解为 x R ∖ ( Q ∖ b ) xR\setminus (Q\setminus b) xR∖(Q∖b) 或 x E ∖ ( R ∖ ( Q ∖ b ) ) 。 xE\setminus (R\setminus (Q\setminus b))。 xE∖(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. ⎩ ⎨ ⎧2x1x2−5x3x413x1−5x27x4−92x2x3−x46x16x2−x3−4x40程序如下 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);xR\(Q\b)程序运行结果如下
x -66.555625.6667-18.777826.5556
或采用 QR 分解的第二种格式程序如下 [Q,R,E]qr(A);xE*(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 XRR XR′R。MATLAB 函数 chol(X) 用于对矩阵 X X X 进行 Cholesky 分解其调用格式如下。① Rchol(X)产生一个上三角阵 R R R使 R ′ R X RRX R′RX。若 X X X 为非对称正定则输出一个出错信息。② [R,p]chol(X)这个命令格式将不输出出错信息。若 X X X 为对称正定的则 p 0 p0 p0 R R R 与上述格式得到的结果相同否则 p p p 为一个正整数。如果 X X X 为满秩矩阵则 R R R 为一个阶数为 q p − 1 qp-1 qp−1 的上三角阵且满足 R ′ R X ( 1 : q , 1 : q ) RR X(1:q,1:q) R′RX(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];Rchol(A)R 1.4142 0.7071 0.70710 1.2247 -1.22470 0 1.0000
可以验证 R ′ R A RR A R′RA输入如下程序 R*Rans 2.0000 1.0000 1.00001.0000 2.0000 -1.00001.0000 -1.0000 3.0000
说明结果是正确的。此时我们利用第二种格式对矩阵 A A A 进行 Cholesky 分解 [R,p]chol(A)R 1.4142 0.7071 0.70710 1.2247 -1.22470 0 1.0000p 0
结果中出现 p 0 p0 p0这表示矩阵 A A A 是一个正定矩阵。如果试图对一个非正定矩阵进行 Cholesky 分解则将得出错误信息所以chol 函数还可以用来判定矩阵是否为正定矩阵。实现 Cholesky 分解后线性方程组 A x b Axb Axb 的解为 R ′ R x b RRxb R′Rxb所以 x R ∖ ( R ′ ∖ b ) 。 xR\setminus (R\setminus b)。 xR∖(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. ⎩ ⎨ ⎧2x1x2−5x3x413x1−5x27x4−92x2x3−x46x16x2−x3−4x40程序如下 A[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b[13,-9,6,0];Rchol(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−x29−x110x2−2x37−2x210x36将方程改写为 { 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}(62x_{2}) \end{matrix}\right. ⎩ ⎨ ⎧x110x2−2x3−7x210x1−9x3101(62x2)这种形式的好处是将一组 x x x 代入右端可以立即得到另一组 x x x。如果两组 x x x 相等那么它就是方程组的解不等时可以继续迭代。例如选取初值 x 1 x 2 x 3 0 x_{1}x_{2}x_{3}0 x1x2x30则经过一次迭代后 得到 x 1 − 7 x_{1}-7 x1−7 x 2 − 9 x_{2}-9 x2−9 x 3 0.6 x_{3}0.6 x30.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}^{(k1)}10x_{2}^{(k)}-2x_{3}^{(k)}-7 \\x_{2}^{(k1)}10x_{1}^{(k)}-9 \\x_{3}^{(k1)}\frac{1}{10}(62x_{2}^{(k)}) \end{matrix}\right. ⎩ ⎨ ⎧x1(k1)10x2(k)−2x3(k)−7x2(k1)10x1(k)−9x3(k1)101(62x2(k))
2.1 Jacobi 迭代法
对于线性方程组 A x b Axb Axb如果 A A A 是非奇异方阵且 a i i ≠ 0 ( i 1 , 2 , … , n ) a_{ii}\ne 0(i1,2,…,n) aii0(i1,2,…,n)则可将 A A A 分解为 A D − L − U AD-L-U AD−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 Axb Axb 转化为 x D − 1 ( L U ) x D − 1 b xD^{-1}(LU)xD^{-1}b xD−1(LU)xD−1b与之对应的迭代公式为 x k 1 D − 1 ( L U ) x k D − 1 b x^{k1}D^{-1}(LU)x^{k}D^{-1}b xk1D−1(LU)xkD−1b这就是 Jacobi 迭代公式。如果序列 { x k 1 } \left \{x^{k1}\right \} {xk1} 收敛于 x x x则 x x x 必是方程 A x b Axb Axb 的解。Jacobi 迭代法的 MATLAB 函数文件 jacobi.m 如下
function [y,n]jacobi(A,b,x0,ep)
if nargin3ep1.0e-6;
elseif nargin3errorreturn
end
Ddiag(diag(A)); %求A的对角矩阵
L-tril(A,-1); %求A的下三角阵
U-triu(A,1); %求A的上三角阵
BD\(LU);
fD\b;
yB*x0f;
n1; %迭代次数
while norm(y-x0)epx0y;yB*x0f;nn1;
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−x29−x110x2−2x37−2x210x36在程序中我们调用函数文件 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.99580.95790.7916n 11
2.2 Gauss-Serdel 迭代法
在 Jacobi 迭代过程中计算 x i ( k 1 ) x^{(k1)}_{i} xi(k1) 时 x 1 ( k 1 ) , … , x i − 1 ( k 1 ) x^{(k1)}_{1},…,x^{(k1)}_{i-1} x1(k1),…,xi−1(k1) 已经得到不必再用 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^{(k1)}(LU)X^{(k)}b Dx(k1)(LU)X(k)b 可以改进为 D x ( k 1 ) L ( k 1 ) X ( k ) b Dx^{(k1)}L^{(k1)}X^{(k)}b Dx(k1)L(k1)X(k)b于是得到 x ( k 1 ) ( D − L ) − 1 U x ( k ) ( D − L ) − 1 b x^{(k1)}(D-L)^{-1}Ux^{(k)}(D-L)^{-1}b x(k1)(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 nargin3ep1.0e-6;
elseif nargin3errorreturn
end
Ddiag(diag(A)); %求A的对角矩阵
L-tril(A,-1); %求A的下三角阵
U-triu(A,1); %求A的上三角阵
G(D-L)\U;
f(D-L)\b;
yG*x0f;
n1; %迭代次数
while norm(y-x0)epx0y;yG*x0f;nn1;
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−x29−x110x2−2x37−2x210x36在程序中我们调用函数文件 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.99580.95790.7916n 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\\ 12 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 -27268n 4 [x,n]gauseidel(a,b,[0;0;0])x NaNNaNNaNn 1012
可见对此方程用 Jacobi 迭代收敛而 Gauss-Serdel 迭代不收敛。因此在使用迭代法时要考虑算法的收敛性。
3. 求线性方程的通解
线性方程组的求解分为两类一类是求方程组的唯一解即特解另一类是求方程组的无穷解即通解。这里对线性方程组 A x b Axb Axb 的求解理论进行归纳。1 当系数矩阵 A A A 是一个满秩方阵时方程 A x b Axb Axb 称为恰定方程方程有唯一 解 x A − 1 b xA^{-1}b xA−1b这是最基本的一种情况。一般用 x A ∖ b xA\setminus b xA∖b 求解速度更快。2 当方程组右端向量 b 0 b0 b0 时方程称为齐次方程组。齐次方程组总有零解因此称解 x 0 x0 x0 为平凡解。当系数矩阵 A A A 的秩小于 n n n n n n 为方程组中未知变量的个数时齐次方程组有无穷多个非平凡解其通解中包含 n n n-rank(4) 个线性无关的解向量用 MATLAB 的函数 null(A,r) 可求得基础解系。3 当方程组右端向量 b ≠ 0 b≠0 b0 时系数矩阵的秩 rank(4) 与其增广矩阵的秩 rank([A,b]) 是判断其是否有解的基本条件。① 当 rank(A)rank([4,b])n 时方程组有唯一解 x r A ∖ b xrA\setminus b xrA∖b 或 x p i n v ( A ) ∗ b xpinv(A)*b xpinv(A)∗b。② 当 rank(4)rank([A,b])n 时方程组有无穷多个解其通解 方程组的一个特解 对应的齐次方程组 A x 0 Ax0 Ax0 的通解。可以用 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[];
if norm(b)0 %非齐次方程组if rank(A)rank([A,b])if rank(A)n %有唯一解disp(原方程组有唯一解x);xA\b; else %方程组有无穷多个解基础解系disp(原方程组有无穷个解特解为x,其齐次方程组的基础解系为y);xA\b; ynull(A,r);endelsedisp(方程组无解); %方程组无解x[];end
else %齐次方程组disp(原方程组有零解x) ;xzeros(n,1); %0解if rank(A)ndisp(方程组有无穷个解基础解系为y); %非0解ynull(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−2x23x3−x413x1−x25x3−3x422x1x22x3−2x43程序如下 A[1,-2,3,-1;3,-1,5,-3;2,1,2,-2];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. ⎩ ⎨ ⎧x1x2−3x3−x413x1−x2−3x34x44x15x2−9x3−8x40程序如下 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,yformat short %恢复默认的短格式输出程序运行结果如下
警告: 秩亏秩 2tol 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 ] Xk_{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} Xk1 3/23/210 k2 −3/47/401 00−8/153/5 其中 k 1 、 k 2 k_{1}、k_{2} k1、k2 为任意常数。
