学习MATLAB线性系统与线性方程

Today:
• Linear equation 线性方程
• Linear system 线性系统

Gaussian Elimination –rref() 高斯消去法

例子

代码

>> A = [1 2 1;2 6 1;1 1 4];
b = [2; 7; 3];
R = rref([A b])

运行结果

LU Factorization LU分解

LU分解步骤

分解后得到的L和U分别是上三角矩阵和下三角矩阵

How to ObtainL and U?

LU Factorization Example

怎样在MATLAB中实现LU分解呢？
还是以上面的矩阵为例

代码

A = [1 1 1;2 3 5;4 6 8];
[L, U, P] = lu(A)

运行结果

解释
[L,U] = lu(A) 将满矩阵或稀疏矩阵 A 分解为一个上三角矩阵 U 和一个经过置换的下三角矩阵 L，使得 A = L*U。
[L,U,P] = lu(A) 还返回一个置换矩阵 P，并满足 A = P’LU。在此语法中，L 是单位下三角矩阵，U 是上三角矩阵。

Matrix Left Division: \or mldivide()矩阵左除法

代码及运行

>>  A=[1 2 1;2 6 1;1 1 4];
>> b=[2;7;3];
>> x=A\b

结果

练习

代码

A=[R1 0 0 R4 0;0 R2 0 -R4 R5;0 0 -R3 0 R5; 1 -1 0 -1 0;0 1 -1 0 -1];
>> 
>> 
>> 
>> syms R1 R2 R3 R4 R5 V1 V2
A=[R1 0 0 R4 0;0 R2 0 -R4 R5;0 0 -R3 0 R5; 1 -1 0 -1 0;0 1 -1 0 -1];
>> b=[V1;0;V2;0;0]
 
b =
 
V1
 0
V2
 0
 0
 
>> x=A\b

运行结果

Matrix Decomposition Functions 矩阵分解函数

Cramer’s (Inverse) Method 克拉默法则

结论是红色部分

Inverse Matrix 逆矩阵

Solving Equations Using Cramer’s Method

用克拉默法则去解方程组

代码

>> A=[1 2 1;2 6 1;1 1 4];
>> b=[2;7;3];
>> x=inv(A)*b

运行结果

但是克拉默法则并不总是起作用，因为可能存在矩阵A不可逆的情况
即奇异矩阵的情况

Singular奇异

以下两种标红的情况在解方程时无解或有无穷多解，此情况被称为奇异
此时不能使用克拉默法则

例子

这两行成比例，不满秩，奇异，行列式为零，逆矩阵不存在
代码

>> A=[1 2 3 4;2 4 6 8;9 8 7 6;1 3 2 8];
>> inv(A)

>> det(A)

运行结果

Problem with Cramer’s Method

Functions to Check Matrix Condition 检查矩阵

检查矩阵代码

>>A = [ 1 2 3; 2 4.0001 6; 9 8 7];
cond(A)

>> B = [ 1 2 3; 2 5 6; 9 8 7];
cond(B)

运行结果

矩阵的逆运算的条件数用于测量线性方程组的解对数据错误的敏感程度。它指示矩阵求逆结果和线性方程解的精度
大的条件数表示系数矩阵 A 中的小小变化会导致线性方程 Ax = b 和 xA = b 的输出 b 发生较大变化。
极端情况是 A 因过于病态而成为奇异矩阵（无限条件数），此时它没有逆矩阵，而线性方程没有唯一解

以上是线性方程的部分，下面来学习线性系统

Linear System 线性系统

Eigenvalues and Eigenvectors 特征值和特征向量

特征值的作用（物理意义）
特征值可以用来放大或缩小特征向量，如果此时输入是用特征向量表示的则，此时特征值和特征向量可以决定输入中的哪一部分被放大，哪一部分又被缩小。这也就是系统矩阵（系统函数）起作用的原因。

举例来理解

当系统矩阵为A时，其特征值为λ1，λ2.对应的特征向量分别是V1,V2.
特征向量V1的方向为图中粉红色箭头向量的方向
特征向量V2的方向为图中粉蓝色箭头向量的方向
所以当输入为上图右侧图片内容时，在蓝色箭头方向上的向量输出是放大三倍，在粉红的箭头方向上的向量输出为不变
对于其他方向上的向量，可以将他们按照V1，V2的方向进行分解，然后再计算输出。

此时的输出为

Solving Eigenvalues and Eigenvectors

以上是数学计算的过程。

由此我们产生一个问题，Eigenvalues and Eigenvectors既然这么好用，怎样用MATLAB求解呢？

eig() 解出特征值和特征向量

代码

>> [v,d]=eig([2 -12;1 -5])

运行结果

Matrix Exponential:expm() 矩阵指数

代码

A =[0 -6 -1; 6 2 -16; -5 20 -10];
x0 = [1 1 1]';
X = [];
for t = 0:.01:1
    X = [X expm(t*A)*x0];
end
plot3(X(1,:),X(2,:),X(3,:),'-o');
xlabel('x_1');
ylabel('x_2');
zlabel('x_3');
grid on;
axis tight square

运行结果

expm 矩阵指数

• Y = expm(X) 计算 X 的矩阵指数。虽然不按此种方式计算，但是如果 X 包含一组完整的特征向量 V 和对应特征值 D，则 [V,D] = eig(X) 且
• expm(X) = V*diag(exp(diag(D)))/V 对于逐个元素的指数运算，请使用 exp。

计算和比较 A 的指数和 A 的矩阵指数。

A = [1 1 0; 0 0 2; 0 0 -1];
exp(A)

结果

expm(A)

结果

两种结果的对角线元素相等，这对于任何三角矩阵都是如此。非对角线元素（包括对角线下方的元素）则不相同。