B站台湾大学郭彦甫｜MATLAB 学习笔记｜12 线性方程式和线性系统 Linear equations

MATLAB学习笔记（12 线性方程式和线性系统 Linear equations）

如想获得更好得阅读体验，可以转到后面链接 12

1. 线性方程 (Linear Equation)

1.1 提出问题

题目：

对于一个电路网络：

给出电压 $V_1$ 和 $V_2$ 、电阻 $R_1...R_5$，求解电流 $i_1...i_5$.

解：

根据基尔霍夫电压电流定律（KCL/KVL），可得到:

[基尔霍夫定律内容](基尔霍夫定律（电学定律）_百度百科 (baidu.com))

$$\begin{array}{rl}& V_1=R_1{i}_1+R_4{i}_4\\ & R_4{i}_4=R_2{i}_2+R_5{i}_5\\ & R_5{i}_5=R_3{i}_3+V_2\\ & i_1=i_2+i_4\\ & i_2=i_3+i_5\end{array}$$
转换为矩阵形式：
$$\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{ccccc}{R}_1& 0& 0& R_4& 0\\ 0& R_2& 0& -R_4& R_5\\ 0& 0& -R_3& 0& R_5\\ 1& -1& 0& -1& 0\\ 0& 1& -1& 0& -1\end{array}\right]}}\underset{\mathbf{x}}{\underbrace{\left[\begin{array}{c}{i}_1\\ i_2\\ i_3\\ i_4\\ i_5\end{array}\right]}}=\underset{\mathbf{b}}{\underbrace{\left[\begin{array}{c}{V}_1\\ 0\\ V_2\\ 0\\ 0\end{array}\right]}}$$

得到的矩阵方程中，$A$ 和 $b$ 已知，$x$ 未知，需要求解 $x$.

1.2 求解线性方程组

可以使用的方法：

1. 消元法 (Successive elimination (通过因式分解))

   包括高斯消去法、LU分解

2. 克莱默法则 (Cramer’s method)

1.3 高斯消去法 (Gaussian Elimination)

1.3.1 数学解释：

将增广矩阵进行 行与行之间的化简，尽可能使增广矩阵的左下方数字为 0

example:

给定方程组：

$$\{\begin{array}{r}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}$$
求解过程：
$$\begin{array}{rl}& \{\begin{array}{r}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}⟹\left[\begin{array}{cccc}1& 2& 1& 2\\ 2& 6& 1& 7\\ 1& 1& 4& 3\end{array}\right]\text{（写成增广矩阵的形式）}\\ & \Rightarrow \left[\begin{array}{cccc}1& 2& 1& 2\\ 0& 2& -1& 3\\ 0& -1& 3& 1\end{array}\right]\text{（行与行相减）}\Rightarrow \left[\begin{array}{cccc}1& 2& 1& 2\\ 0& 2& -1& 3\\ 0& 0& 5/2& 5/2\end{array}\right]\end{array}$$
​ 通过最后的矩阵列方程组求出解 $x=[x_1,x_2,x_3]$

1.3.2 MTALAB 方法：rref()

A = [1 2 1;2 6 1;1 1 4];
b = [2; 7; 3];
R = rref([A b])	%使用 rref 以简化行阶梯形矩阵表示该方程组。

>> Class12

R =

     1     0     0    -3
     0     1     0     2
     0     0     1     1

1.4 LU分解 (LU Factorization)

• 假设我们要求解： $Ax=b$, 其中 $A\in R^{m\times m}$

• 将 $A$ 分解为2个三角矩阵: $A=L^{-1}U$

• 则问题转化为：$Ax=b⟹L^{-1}\underset{y}{\underbrace{Ux}}=b$

• 策略方法：

1. 求解 $L^{-1}y=b$ 得到 $y$
2. 求解 $Ux=y$

1.4.1 上/下三角矩阵定义

• 下三角矩阵 $\mathbf{L}=\left[\begin{array}{ccc}1& 0& 0\\ ⋮& \ddots & 0\\ ⋮& \cdots & 1\end{array}\right]\in {\mathfrak{R}}^{m\times m}$（对角线为1）

• 上三角矩阵 $\mathbf{U}=\left[\begin{array}{ccc}⋮& \cdots & ⋮\\ 0& \ddots & ⋮\\ 0& 0& ⋮\end{array}\right]\in {\mathfrak{R}}^{m\times m}$

1.4.2 如何获得 $L$ 和 $U$ ?

​ 使用一系列的左乘：$
\underbrace{\boldsymbol{L}{m} \ldots \boldsymbol{L}{2} \boldsymbol{L}{1}}{\boldsymbol{L}} \boldsymbol{A}=\boldsymbol{U}$

1.4.2 LU 分解举例

对于一个矩阵：$\mathbf{A}=\left[\begin{array}{lll}1& 1& 1\\ 2& 3& 5\\ 4& 6& 8\end{array}\right]$

解：

​ 先左乘一个矩阵 $L_1$，其中 $L_1$ 是在单元矩阵 $\mathbf{E}=\left|\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right|$ 的基础上补充得到，$L_{1}A$ 相乘使得矩阵 $A$ 中的元素 $x_{21}=2$ 和 $x_{31}=4$ 化简为 $0$：
$${\mathbf{L}}_1\mathbf{A}=\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -4& 0& 1\end{array}\right]\left[\begin{array}{lll}1& 1& 1\\ 2& 3& 5\\ 4& 6& 8\end{array}\right]=\left[\begin{array}{lll}1& 1& 1\\ 0& 1& 3\\ 0& 2& 4\end{array}\right]$$
​ 同理，左乘一个矩阵 $L_2$ ，$L_2(L_{1}A)$ 使得矩阵 $L_{1}A$ 中的元素 $x_{32}=2$ 化简为 $0$。
$${\mathbf{L}}_2\left({\mathbf{L}}_1\mathbf{A}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& -2& 1\end{array}\right]\left[\begin{array}{lll}1& 1& 1\\ 0& 1& 3\\ 0& 2& 4\end{array}\right]=\left[\begin{array}{rrr}1& 1& 1\\ 0& 1& 3\\ 0& 0& -2\end{array}\right]=\mathbf{U}$$

由此求解出 $L=L_2{L}_1$ 和 $U$.

1.4.3 MATLAB 方法: lu()

• matlab 中的 lu() 函数在输入一个矩阵 A 之后,可以得到 L, U, P，其语法为 [L,U,P]=lu(A) ，数学关系满足：$PA=LU$
• 计算 $Ax=b$, 由 $PA=LU$ 可得到: $P^{-1}LUx=b$, $\text { 令 } y=L^{-1} Pb \text {, 有 } U x=y $

example:

给出：
$$\mathbf{A}=\left[\begin{array}{lll}1& 1& 1\\ 2& 3& 5\\ 4& 6& 8\end{array}\right]\mathbf{b}=\left[\begin{array}{l}2\\ 7\\ 3\end{array}\right]$$

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

• 方法一：(matlab 中的数学计算方法)

​ 得到 $L,P,b$ 可利用 $ y=L^{-1} Pb$ 求 $y$，可利用 $Ux=yb$ 关系式求解 $x$。

• 方法二：(老师的方法)
  $$\{\begin{array}{c}{L}^{-1}y=b\\ Ux=y\end{array}$$

1.4.4 矩阵的逆的快速表示：\

对于方程组：
$$\{\begin{array}{r}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}$$

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

exercise (P14)（有疑惑）

题目：

$\text{ Write a function to solve }{i}_1\dots i_5\text{ for given }{V}_1,V_2\text{, and }{R}_1\dots R_5$
$$\begin{array}{rl}& V_1=R_1{i}_1+R_4{i}_4\\ & R_4{i}_4=R_2{i}_2+R_5{i}_5\\ & R_5{i}_5=R_3{i}_3+V_2\\ & i_1=i_2+i_4\\ & i_2=i_3+i_5\end{array}$$

解：

% LU分解
clear
clc
syms R1 R2 R3 R4 R5 i1 i2 i3 i4 i5 V1 V2 V3 V4 V5
%根据方程组列出矩阵方程Ax=b,其中x与电流i有关
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];

x=[i1; i2; i3; i4; i5];
b=[V1; V2; V3; V4; V5];
[L,U,P]=lu(A);

% matlab的数学计算方法
x1 = x;
y1 = L\(P*b);
x1 = U\y1；


% 老师的方法
x2 = x;
y2 = L*b;
x2 = U\y2；

（解出来的 x1 和 x2 不一致，比较疑惑）

1.5 克莱默法则（Cramer’s (Inverse) Method）

内容：

• 问题如下:

$$\underset{A}{\underbrace{\left[\begin{array}{cc}3& -2\\ 1& 4\end{array}\right]}}\underset{x}{\underbrace{\begin{array}{l}x\\ y\end{array}]}}=\underset{b}{\underbrace{\left[\begin{array}{c}5\\ 11\end{array}\right]}}$$

• 假设存在 $A^{-1}$ 使得 $A{A}^{-1}=A^{-1}A=E$
• 则变量 $x$ 为：$x=A^{-1}b$

1.5.1 逆矩阵（Inverse Matrix）

• 对于矩阵 $A$，其逆矩阵定义如下：

$${\mathbf{A}}^{-1}={\left[\begin{array}{ll}a& b\\ c& d\end{array}\right]}^{-1}=\frac{1}{\det (\mathbf{A})}\mathrm{adj}(\mathbf{A})=\frac{1}{\det (\mathbf{A})}\left[\begin{array}{cc}d& -b\\ -c& a\end{array}\right]$$

​ 其中 $det(A)$ 表示 $A$ 的行列式：$\det (\mathbf{A})=|ad-bc|$

• 性质：$\mathbf{A}=({\mathbf{A}}^{-1}{)}^{-1},(k\mathbf{A}{)}^{-1}=k^{-1}{\mathbf{A}}^{-1}$

example:（使用克莱姆法则求解方程组）

题目：
$$\{\begin{array}{r}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}\Rightarrow \left[\begin{array}{lll}1& 2& 1\\ 2& 6& 1\\ 1& 1& 4\end{array}\right]x=\left[\begin{array}{l}2\\ 7\\ 3\end{array}\right]$$

解：

$x=A^{-1}b$

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

exercise (P19)

题目：

• 在三维空间画该平面 $\{\begin{array}{c}x+y+z=0\\ x-y+z=0\\ x+3z=0\end{array}$

解：

clear
clc
clf
x = 0:40;
y = 0:40；
[X1,Y1]=meshgrid(x,y);	%画三维图之前绘制网格，可以参考高阶绘图部分讲义
[X2,Y2]=meshgrid(x,y);
[X3,Y3]=meshgrid(x,y);

Z1 = -X1-Y1;
Z2 = -X2+Y2;
Z3 = -1/3*X3;

hold on
surf(Y1,X1,Z1);
surf(Y2,X2,Z2);
surf(Y3,X3,Z3);
axis square;
hold off

%如果需要改变Y轴的方向和z轴的取值范围，可以在 Figure-查看-属性检查器 中更改

1.6 检查矩阵状况的函数：cond()

如果矩阵 $A$ 出现下图中第二和第三种情况，则 $A$ 的行列式为零（$det(A)=0$，Singular）：

黑色文字不用理，是截图问题

在 MATLAB 中可以使用 cond() 函数观察矩阵的“健康状况”，得到的值越小，矩阵对求逆运算越不敏感。

2. 线性系统 (Linear System)

2.1 提出问题

线性方程式是 $Ax=b$ 的形式，线性系统是 $Ab=y$ 的形式

• 给出方程组：$\{\begin{array}{l}2\cdot 2-12\cdot 4=x\\ 1\cdot 2-5\cdot 4=y\end{array}$

• 矩阵表达形式：

2.2 特征值和特征向量 (Eigenvalues and Eigenvectors)

2.2.1 数学方法：

• 对于一个系统 $A\in R^{m\times m}$，如果矩阵相乘 $y=Ab$ 比较复杂，可以找到一个向量 $v$，使得 $\mathbf{A}{\mathbf{v}}_i={\lambda }_i{\mathbf{v}}_i$ ，并且 ${\lambda }_i$ 为常数（把 $v_i$ 称为特征向量，${\lambda }_i$ 称为特征值）

• 然后分解 $\mathbf{b}=\sum {\alpha }_i{\mathbf{v}}_i,{\alpha }_i\in \Re$

• 矩阵相乘变为：

$$\mathbf{A}\mathbf{b}=\sum {\alpha }_i\mathbf{A}{\mathbf{v}}_i=\sum {\alpha }_i{\lambda }_i{\mathbf{v}}_i$$

2.2.2 特征值和特征向量的解释

对于矩阵 $A$ 的特征值和特征向量：
$$\{\begin{array}{l}A{v}_1={\lambda }_1{v}_1\\ A{v}_2={\lambda }_2{v}_2\end{array}$$
​ 当输入为 $v_1$ 时, 如果 ${\lambda }_1>1$, 称系统对输入进行放大；如果输入不是 $v_1$, 可以对输入进行分解, 输入 $=\left({\alpha }_1{v}_1+{\alpha }_2{v}_2\right)$

example:

• 给出 $\mathbf{A}\mathbf{b}=\left[\begin{array}{cc}2& -12\\ 1& -5\end{array}\right]\left[\begin{array}{l}2\\ 4\end{array}\right]$ （假如不想计算矩阵乘法，用求特征值的方法）

• 利用 MATLAB 的 eig() 函数可以得到下面结果：（2.2.3部分有解释）

$${\lambda }_1{\mathbf{v}}_1=\left[\begin{array}{l}0.97\\ 0.24\end{array}\right],{\lambda }_2{\mathbf{v}}_2=-2\left[\begin{array}{l}0.95\\ 0.32\end{array}\right]$$

• 把 $b$ 拆解 $\mathbf{b}={\alpha }_1{\mathbf{v}}_1+{\alpha }_2{\mathbf{v}}_2=-41.2{\mathbf{v}}_1+44.3{\mathbf{v}}_2$

• 计算 $Ab$

$$\begin{gathered} \mathbf{Ab}=\mathbf{A}\left({\alpha }_1{\mathbf{v}}_1+{\alpha }_2{\mathbf{v}}_2\right) \\ ={\alpha }_1{\mathbf{Av}}_1+{\alpha }_2\mathbf{A}({\alpha }_1={\mathbf{v}}_1+{\alpha }_2{\mathbf{v}}_1{\mathbf{v}}_1+{\alpha }_2{\lambda }_2{\mathbf{v}}_2 \\ =(-41.2)(-1)\left[\begin{array}{l}0.97\\ 0.24\end{array}\right]+(44.3)(-2)\left[\begin{array}{l}0.95\\ 0.32\end{array}\right] \end{gathered}$$

2.2.3 算特征值和特征向量用 eig() 函数

求特征值和特征向量：
$$A=\left[\begin{array}{cc}2& -12\\ 1& -5\end{array}\right]$$

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

2.3 矩阵指数：expm()