matlab语言与应用 04 线性代数

现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

04 线性代数问题的计算机求解

04.01 特殊矩阵输入

零矩阵、幺矩阵及单位矩阵
A = zeros(n),
B = ones(n),
C = eye(n)
A = zeros(m, n)
B = ones(m, n)
C = eye(m, n)
A = zeros(size(B))
A = zeros(n, m, …, k)

例4-1 特殊矩阵的建立

% 生成一个3x8的零矩阵A，生成一个和A维数相同的单位矩阵B
A = zeros(3, 8),
B = eye(size(A))
% 注意: zeros() 和 ones() 也可用于多为数组的生成

随机元素矩阵
均匀分布随机数矩阵的输入
[0, 1]区间，A = rand(n), A = rand(n, m)
[a, b]区间, B = a + (b-a)*rand(n, m)
征态分布随机数矩阵的输入
标准征态分布N(0, 1)： A = randn(n), A = randn(n, m)
$N(\mu, \sigma^2): B = \mu + \sigma* \mathbb{randn}(n, m)$

对角元素矩阵
对角矩阵的数学描述
$\mathbb{diag}(a_1, a_2, \cdots, a_n) = \begin{bmatrix}a_1 & & & \\ & a_2 & & \\ & & \ddots & \\ & & & a_n \end{bmatrix}$
其中,所有的非对角元素都为0

diag()函数的使用
已知向量生成对角元素: A = diag(V)
已知矩阵提取对角元素列向量: V = diag(A)
生成主对角线上第k条对角线为v的矩阵: A = diag(v, k)
注意: k可为负整数

例4-2 diag()函数
diag()函数的不同调用格式

% 生成对角矩阵
C = [1 2 3];
V = diag(C)
% 对角元素提取
V1 = diag(V)
% 生成主对角线上第2条对角线为C的矩阵
C = [1 2 3];
V2 = diag(C, 2), 
V3 = diag(C, -1)

Hankel矩阵
Hankel矩阵的一般形式
$H = \begin{bmatrix}c_1 & c_2 & \cdots & c_m \\ c_2 & c_3 & \cdots & c_{m+1} \\ \vdots & \vdots & \ddots & \vdots \\ c_n & c_{n+1} & \cdots & c_{n+m-1} \end{bmatrix}$
反对角线上元素相同
函数调用:
H = hankel(C)
H = hankel(C, R)
C(end) = R(1)

Vandermonde矩阵
Vandermonde矩阵的数学描述
$V = \begin{bmatrix}c_1^{n-1} & c_2^{n-2} & \cdots & c_1 & 1 \\ c_2^{n-1} & c_2^{n-2} & \cdots & c_2 & 1\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ c_n^{n-1} & c_n^{n-2} & \cdots & c_n & 1 \end{bmatrix}$
其中，$v_{i,j} = c_i^{n-j}, i,j = 1, 2, \cdots, n$
生成Vandermonde矩阵 V = vander(c)

相伴矩阵
一个首一化的多项式
$P(s) = s^n + a_1s^{n-1} + a_2s^{n-2} + \cdots + a_{n-1}s + a_n$
其相伴矩阵
$A_c = \begin{bmatrix}-a_1 & -a_2 & \cdots & -a_{n-1} & -a_n \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix}$
matlab函数: B = compan(p)

例4-5 由多项式构造相伴矩阵
考虑一个多项式 $P(s) = 2s^4 + 4s^2 + 5s + 6$
求其相伴矩阵
非首一多项式，可以自动转换

P = [2 0 4 5 6];
A = compan(P)

符号矩阵输入
数值矩阵A转符号矩阵B: B = sym(A)
早期版本支持重载函数的方式，并置于目录@sym下，可编写compan, hankel, vander
最新版本这三个函数直接支持符号矩阵

任意矩阵输入
任意向量输入方法
行向量: A = sym(‘a%d’, [1, n])
列向量: A = sym(‘a%d’, [n, 1])
任意矩阵的输入方法
方阵输入: A = sym(‘a%d%d’, n)
任意矩阵输入: A = sym(‘a%d%d’, [n, m])

例4-7 任意矩阵的生成
输入如下任意矩阵
$A = \begin{bmatrix}a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}$
$B = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32} \\ a_{41} & a_{42} \end{bmatrix}$
$C = \begin{bmatrix}f_{11} & f_{12} & f_{13} & f_{14} \\ f_{21} & f_{22} & f_{23} & f_{24} \\ f_{31} & f_{32} & f_{33} & f_{34} \\ f_{41} & f_{42} & f_{43} & f_{44} \end{bmatrix}$

A = sym('a%d%d', 4),
B = sym('a%d%d', [4,2]),
C = sym('f%d%d', [4, 4])

稀疏矩阵输入
稀疏矩阵:矩阵大部分元素为零,仅少部分元素非零.
稀疏矩阵输入: A = sparse(p, q, w)
稀疏矩阵与常规矩阵转换：B = full(A), A = sparse(B)

04.02 矩阵性质分析

行列式
矩阵$A = \lbrace a_{ij} \rbrace$ 的行列式定义为
$D = |A| = det(A) = \sum (-1)^k a_{1k_1}a_{2k_2}\cdots a_{nk_n}$
函数调用:
d = det(A), d = det(sym(A))
注意: 该函数可用于数值和符号运算，方法和A的类型一致
自动识别矩阵的形式，选择求解方法

例4-8 矩阵求行列式
$A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{bmatrix}$

A = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
det(A)
det(sym(A))

例4-10 任意矩阵的行列式
任意4x4矩阵的行列式

A = sym('a%d%d', 4);
d = det(A)

代数余子式$A_{23}^*$

i = 2;
j = 3;
B = A;
B(i, :) = [];
B(:, j) = [];
A23 = (-1)^(i+j)*det(B)
% 下面语句也可求得代数余子式
syms a23;
A23_1 = simplify((d-subs(d, a23, 0))/a23)

矩阵的迹
假设一个方阵为 $A = \lbrace a_{ij} \rbrace, i,j = 1, 2, \cdots, n$
矩阵的迹定义为对角线元素之和：tr(A) = $\sum \limits_{i=1}^{n}a_{ij}$
函数调用格式：t = trace(A)
迹也等于特征值的和：sum(diag(A))

矩阵的秩
矩阵A的秩:线性不相关的最大行(列)数
rank(A) $= r_c = r_r$
数值方法或符号方法: r = rank(A)
给定精度$\varepsilon$下求数值秩: r = rank(A, $\varepsilon$)

向量的范数
范数式测度 – 用一个值描述向量、矩阵的大小
$函数\rho(x)为x向量的范数的条件:$
$\rho(x) \geq 0 且 \rho(x) = 0 的充要条件是x = 0$
$\rho(ax) = |a| \rho(x),a 为任意标量$
$对向量x和y，有\rho(x+y) \leq \rho(x) + \rho(y)$
向量的p-范数
$\Vert x \Vert _{\infty} = \max \limits_{1 \leq i \leq n} |x_i|$
$\Vert x \Vert _{p} = (\sum \limits_{i = 1}^{n} |x_i|^p)^{\frac{1}{p}}, p = 1, 2, \cdots,$

矩阵的范数
对于任意非零向量x, 矩阵A的范数是
$\Vert A \Vert = \sup \limits_{x \neq 0} \dfrac{\Vert A_x \Vert}{\Vert x \Vert}$
其它常用范数
$\Vert A \Vert _1 = \max \limits_{1 \leq j \leq n} \sum \limits_{i=1}^{n} |a_{ij}|, \Vert A \Vert_2 = \sqrt{s_{max}(A^TA)}$
$\Vert A \Vert _{\infty} = \max \limits_{1 \leq i \leq n} \sum \limits_{i=1}^{n} |a_{ij}|, \Vert A \Vert _F = trace(A^TA)$

矩阵范数计算
s(X)为x矩阵的特征值， $s_{max}(A^TA)$
A矩阵的范数为矩阵AA^T其最大特征值
函数调用格式: $\Vert A \Vert_2 \qquad$ N = norm(A)
选项可以是1，2， inf, ‘fro’: N = norm(A, opts)
注意: 早期版本norm()函数仅用于数值矩阵

矩阵特征多项式
由矩阵A的如下多项式:
$c(s) = det(sI - A) = s^n + c_1s^{n-1} + \cdots + c_{n-1}s + c_n$
多项式c(s)是矩阵A的特征多项式
matlab下多项式由降幂排列的系数向量表示
c = poly(A)
c = charpoly(sym(A))
c = charpoly(sym(A), x)

矩阵多项式求解
矩阵多项式的数学表示:
$B = a_1A^n + a_2A^{n-1} + \cdots + a_nA + a_{n+1}I$
函数调用格式，A不能为符号矩阵
B = polyvalm(a, A)
其中，a = [$a_1, a_2, \cdots, a_n, a_{n+1}$]是特征多项式的按降幂排列的系数

符号运算多项式求值
调用格式: B = polyvalmsym(p, A)

% 符号运算多项式求值
function B = polyvalmsym(p, A)
  E = eye(size(A));
  B = zeros(size(A));
  n = length(A);
  for i = n+1:-1:1
    B = B + p(i)*E;
    E = E*A;
  end;
end

p为多项式系数向量

矩阵多项式点运算
点运算方式定义多项式运算
C = $a_1$x .^n + $a_2$x.^(n-1) + … + $a_{n+1}$
调用格式：C = polyval(a, x)
给出多项式p(符号运算工具箱)：C = subs(p, s, x)

Cayley-Hamilton定理
若矩阵A的特征多项式为
$f(x) = det(sI -A) = a_1s^n + a_2s^{n-1} + \cdots + a_ns + a_{n+1}$
则$f(A) = 0$, 即
$a_1A^n + a_2A^{n-1} + \cdots + a_nA + a_{n+1}I = 0$
可以通过matlab验证任意低阶矩阵

例4-18 验证任意矩阵满足Cayley-Hamilton定理
验证一般5x5矩阵满足Cayley-Hamilton定理

A = sym('a%d%d', 5);
p = charpoly(A);
E = polyvalmsym(p, A);
norm(simplify(E))

04.03 逆矩阵与广义矩阵

矩阵的逆
逆矩阵数学描述 AC = CA = I
其中，A 是n x n 的非奇异方阵，那么 $C = A^{-1}$
函数表示: C = inv(A)
基本行变换(行阶梯型矩阵): $H_1 = \mathbb{rref}(H)$

例4-20 Hilbert矩阵的逆

% 求Hilbert矩阵的逆。
% 4 x4 Hilbert 矩阵数值解
H = hilb(4);
H1 = inv(H),
norm(H*H1 - eye(4))
% 13 x 13 Hilbert 矩阵
H = hilb(13);
H1 = inv(H);
n1 = norm(H*H1 - eye(size(H)))
H2 = invhilb(13);
n2 = norm(H*H2 - eye(size(H)))
% 警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND =  1.334996e-18。

符号运算

% 符号矩阵
% 7 x 7 Hilbert 矩阵
H = sym(hilb(7));
H1 = inv(H)
% 50 x 50 Hilbert 逆矩阵
H = sym(hilb(50));
norm(H*inv(H) - eye(size(H)))

例4-21 奇异矩阵的逆
求如下奇异矩阵的逆
$A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{bmatrix}$

A = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
B = inv(A), 
A*B

使用符号运算工具箱

A = sym(A);
inv(A)
% ans = FAIL

例4-22 符号矩阵的逆
推导如下 Hankel 矩阵的逆
$A = \begin{bmatrix}a_1 & a_2 & a_3 & a_4 \\ a_2 & a_3 & a_4 & 0 \\ a_3 & a_4 & 0 & 0 \\ a_4 & 0 & 0 & 0 \end{bmatrix}$

a = sym('a%d', [1, 4]);
H = hankel(a);
inv(H)

矩阵的广义逆
适用于奇异或矩形矩阵
如果存在ANA = A,则称 N为A矩阵的广义逆矩阵，记作: $N = A^-$
不唯一
定义范数最小化指标
$\min \limits_{N} \Vert AN - I \Vert$

Moore-Penrose广义逆
矩阵M为矩阵A的Moore-Penrose广义逆矩阵的条件:
(i) AMA = A
(ii) MAM = M
(iii) AM 和 MA 均为 Hermite 对称矩阵
记为 $M = A^+$, 又称 伪逆
Moore-Penrose 广义逆矩阵是唯一的
M = pinv(A), M = pinv(A, $\epsilon$)

例4-24 奇异矩阵的伪逆
使用pinv()函数计算矩阵A的伪逆
$A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{bmatrix}$

A = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
B = pinv(A),
A*B

伪逆检验

% 检验 Moore-Penrose 广义逆的三个条件
norm(A*B*A - A),
norm(B*A*B - B),
norm(A*B - (A*B)'),
norm(B*A - (B*A)')
% 检验 (A^+)^+ = A
pinv(pinv(A)), norm(ans - A)

长方形奇异矩阵的伪逆
$A = \begin{bmatrix}6 & 1 & 4 & 2 & 1 \\ 3 & 0 & 1 & 4 & 2 \\ -3 & -2 & -5 & 8 & 4 \end{bmatrix}$

% 求秩并求伪逆
A = [6, 1, 4, 2, 1; 3, 0, 1, 4, 2; -3, -2, -5, 8, 4];
rank(A)
iA = pinv(A), 
norm(A*iA*A - A),
norm(iA*A - A'*iA')
norm(iA*A - A'*iA'),
norm(A*iA - iA'*A')

04.04 特征值与特征向量

一般矩阵的特征值与特征向量
数学描述 $Ax = \lambda x$
非零向量x是特征向量，数值$\lambda$ 是特征值
调用格式: d = eig(A) 或 [V, D] = eig(A)

例4-26 矩阵的特征值
计算下述矩阵的特征值与特征向量
$A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{bmatrix}$

A = [16, 2, 3, 13; 5, 11, 10, 8; 9, 7, 6, 12; 4, 14, 15, 1];
eig(A),
[v, d] = eig(A)
eig(sym(A)),
vpa(ans, 70),
[v, d] = eig(sym(A))

特征向量矩阵的应用
如果矩阵没有重复特征根，则 $V^{-1}AV = D$
特征向量矩阵可以将一般矩阵变换为对角矩阵

A = [16, 2, 3, 13; 5, 11, 10, 8; 9, 7, 6, 12; 4, 14, 15, 1];
[v, d] = eig(sym(A))
s = inv(v)*A*v;
simplify(s)

特征向量不唯一
数值eig函数得出的特征向量是归一化的特征向量

符号矩阵特征值计算
假设已知向量 $c = [1, c_1, c_2, c_3]$
可以构造出相伴矩阵
$H = \begin{bmatrix}-c_1 & -c_2 & -c_3 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}$
试求其特征值与特征向量矩阵

syms c1 c2 c3 real;
H = compan([1 c1 c2 c3]);
[V, D] = eig(H)

矩阵的广义特征向量问题
数学描述: $Ax = \lambda B x$
非零向量x是特征向量, 数值 $\lambda$ 是特征值，
B 是正定对称矩阵
调用格式 d = eig(A, B)
或 [V, D] = eig(A, B)

例4-27 广义特征值
已知A, B,求广义特征值与特征向量矩阵
$A = \begin{bmatrix}5 & 7 & 6 & 5 \\7 & 10 & 8 & 7 \\6 & 8 & 10 & 9 \\5 & 7 & 9 & 10 \end{bmatrix}, B = \begin{bmatrix}2 & 6 & -1& -2 \\5 & -1 & 2 & 3 \\ -3 & -4 & 1 & 10 \\5 & -2 & -3 & 8 \end{bmatrix}$

B = [2, 6, -1, -2; 5, -1, 2, 3; -3, -4, 1, 10; 5, -2, -3, 8];
A = [5, 7, 6, 5; 7, 10, 8, 7; 6, 8, 10, 9; 5, 7, 9, 10];
[V, D] = eig(A, B),
norm(A*V - B*V*D)

04.05 矩阵相似变换与三角分解

矩阵相似变换与正交矩阵
矩阵的相似性的数学描述
相似变换不改变矩阵的特征值
$X = B^{-1}AB$
正交矩阵的数学描述
$Q^{H}Q = I, QQ^{H} = I$
函数调用格式: Q = orth(A)

矩阵的三角分解
一般矩阵分解成下三角与上三角矩阵的积
数学描述: A = LU
其中
$L = \begin{bmatrix}1 & & & \\ l_{21} & 1 & & \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1 \end{bmatrix}, U = \begin{bmatrix}u_{11} & u_{12} & \cdots & u_{1n} \\ & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} \end{bmatrix}$

三角分解的递推算法
递推初值
$u_{1i} = a_{1i}, i = 1, 2, \cdots, n$
递推计算公式
$l_{ij} = \dfrac{a_{ij} - \sum \limits_{k=1}^{j-1} l_{ik}u_{kj}}{u_{jj}}, (j < i)$
$u_{ij} = a_{ij} - \sum \limits_{k=1}^{i-1} l_{ik}u_{kj}, (j \geq i)$

三角分解的matlab求解
lu分解:
数学描述:A = LU
matlab求解: [L, U] = lu(A)
P为置换矩阵
数学表示: $A = P^{-1}LU$
matlab求解: [L, U, P] = lu(A)

例4-30 矩阵的三角分解
给定 $A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{bmatrix}$
两种方法调用lu()函数

A = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
[L1, U1] = lu(A)
% 变换矩阵
[L2, U2, P2] = lu(A)

处理带有主元素的分解

A = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
% 有主元素和置换矩阵的
[L, U, P] = lu(A)
% 验证三角分解公式
inv(P)*L*U
% 三角分解的解析解
A = sym(A);
[L2, U2] = lu(A)

对称矩阵的三角矩阵–Cholesky分解
数学描述: $A = LL^T$
下三角矩阵:
$L = \begin{bmatrix}l_{11} & & & \\ l_{21} & l_{22} & & \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{bmatrix}$
其中 A 是对称矩阵

Cholesky分解的算法
对称矩阵的Cholesky分解算法
$l_{ii} = \sqrt{a_{ii} - \sum \limits_{k=1}^{i-1} l_{ik}^2}$
$l_{ji} = \dfrac{1}{l_{jj}}(a_{ij} - \sum \limits_{k=1}^{j-1} l_{ik}l_{jk}), j < i$
调用格式: D = chol(A)

例4-26 Cholesky分解
对如下A矩阵进行Cholesky分解
$A = \begin{bmatrix}9 & 3 & 4 &2 \\ 3 & 6 & 0 & 7 \\ 4 & 0 & 6 & 0 \\ 2 & 7 & 0 & 9 \end{bmatrix}$

A = [9, 3, 4, 2; 3, 6, 0, 7; 4, 0, 6, 0; 2, 7, 0, 9];
D = chol(A),
D1 = chol(sym(A))

04.06 矩阵Jordan变换与奇异值分解

一般矩阵变成相伴矩阵
如果找到一个列向量x
使得T为非奇异矩阵
$T = [x, Ax, \cdots, A^{n-1}x]$
矩阵A可以转换成一个与类似相伴矩阵的形式
这样的x向量有无穷多个

例4-33 寻找变换矩阵
$A = \begin{bmatrix}5 & 7 & 6 & 5 \\ 7 & 10 & 8 & 7 \\ 6 & 8 & 10 & 9 \\ 5 & 7 & 9 &10 \end{bmatrix}$

A = [5, 7, 6, 5; 7, 10, 8, 7; 6, 8, 10, 9; 5, 7, 9, 10];
while (1)
  x = floor(2*rand(4, 1));
  T = sym([x A*x A^2*x A^3*x]);
  if rank(T) == 4
    break;
  end;
end;
T, A1 = inv(T)*A*T

矩阵对角化
矩阵A特征值互异，则特征向量矩阵T为非奇异矩阵，可将元矩阵变换成对角矩阵
含有复数特征根的矩阵能得出复数的对角矩阵和复数相似变换矩阵
如果含由共轭复数的特征值，则可以将其转换成带有复数块的实数矩阵
如果有复数特征值，则不能将其转换成对角矩阵，只能转换为Jordan矩阵

例4-34 矩阵对角化
试求出下述矩阵的对角矩阵及变换矩阵
$A = \begin{bmatrix}3 & 2 & 2 & 2 \\ 1 & 2 & -2 & -2 \\ -1 & -2 & 0 & -2 \\ 0 & 1 & 3 & 5 \end{bmatrix}$

A = [3, 2, 2, 2; 1, 2, -2, -2; -1, -2, 0, -2; 0, 1, 3, 5];
[v, d] = eig(sym(A));
A1 = inv(v)*A*v

矩阵Jordan变换
Jordan变换用于处理含由重特征值的矩阵
函数调用格式
只返回Jordan矩阵J： J = jordan(A)
返回Jordan矩阵J，和广义特征向量矩阵V
[V, J] = jordan(A)
如果矩阵的特征值互不相同，则函数返回的结果与eig()函数一致

例4-37 Jordan变换
求如下矩阵的Jordan分解
$A = \begin{bmatrix}-71 & -65 & -81 & -46 \\ 75 & 89 & 117 & 50 \\ 0 & 4 & 8 & 4 \\ -67 & -121 & -173 & -58 \end{bmatrix}$

A = [-71, -65, -81, -46; 75, 89, 117, 50; 0, 4, 8, 4; -67, -121, -173, -58];
[V, J] = eig(sym(A))
[V1, J1] = jordan(sym(A))

矩阵的奇异值分解
数学描述: $A^TA \geq 0, AA^T \geq 0$
理论上， $rank(A^TA) = rank(AA^T) = rank(A)$
奇异值定义: $\sigma_i(A) = \sqrt{\lambda_i(A^TA)}$
其中，$\sigma_i$ 是非负特征值

例4-40 奇异值实例
$A = \begin{bmatrix}1 & 1 \\ \mu & 0 \\ 0 & \mu \end{bmatrix} \qquad \mu = 10^{-15}$
求矩阵的秩
$A^TA = \begin{bmatrix}1+\mu^2 & 1 \\ 1 & 1 + \mu^2 \end{bmatrix}$

A = [1, 1; 5*eps, 0; 0, 5*eps];
rank(A)

奇异值分解
奇异值分解 Singular value decomposition(SVD)
$A = LA_1M$
A – n x m 矩阵
L 和 M 为正交矩阵
分解出对角矩阵 $A_1 = diag(\sigma_1, \cdots, \sigma_n)$
对角元素满足: $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_n \geq 0$
求解: S = svd(A)
[L, $A_1$, M] = svd(A)

矩阵的条件数
矩阵的条件数
$cond(A) = \dfrac{\sigma_{\max}}{\sigma_{\min}} = \dfrac{\bar \sigma(A)}{\underline \sigma(A)}$
matlab求解: c = cond(A)
条件数过大，说明矩阵接近奇异值
坏条件矩阵
处理时应该慎重(建议用符号运算)

例4-40 奇异值分解
对下列矩阵进行奇异值分解
$A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{bmatrix}$

A = [16, 2, 3, 13; 5, 11, 10, 8; 9, 7, 6, 12; 4, 14, 15, 1];
[L, A1, M] = svd(A)
% 条件数
cond(A)

例4-41 长方形矩阵分解
对下面的矩阵进行奇异值分解，并验证结果
$A = \begin{bmatrix}1 & 3 & 5 & 7 \\ 2 & 4 & 6 & 8 \end{bmatrix}$

A = [1, 3, 5, 7; 2, 4, 6, 8];
[L, A1, M] = svd(A)
A2 = L*A1 * M',
norm(A-A2)

04.07 线性方程求解

线性方程组求解
数学描述: Ax = B
其中，A和B是给定矩阵
$A = \begin{bmatrix}a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix}, B = \begin{bmatrix}b_{11} & b_{12} & \cdots & b_{1p} \\ b_{21} & b_{22} & \cdots & b_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ b_{m1} & b_{m2} & \cdots & b_{mp} \end{bmatrix}$
线性代数课程：解的三种情况

I.唯一解
唯一解存在条件
当 m = n,且 rank(A) = n, 则唯一解存在
$X = A^{-1}B$
函数调用格式：
X = inv(A)*B
X = A \ B
注意: 推荐使用符号运算方法

例4-42 解方程
求解线性方程组
$\begin{bmatrix}1 & 2 & 3 & 4 \\ 4 & 3 & 2 &1 \\ 1 & 3 &2 & 4\\ 4 &1 & 3 &2 \end{bmatrix} X = \begin{bmatrix}5 & 1 \\ 4 & 2\\ 3 & 3 \\ 2 & 4 \end{bmatrix}$
求解与检验

A = [1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2];
B = [5 1; 4 2; 3 3; 2 4];
x = inv(A) * B, % 数值解
e1 = norm(A*x - B),
x1 = inv(sym(A))*B, % 解析解
e2 = norm(A*x1 - B)

II. 方程组有无穷多解
解的判定矩阵
$C = \begin{bmatrix}a_{11} & a_{12} & \cdots & a_{1n} & b_{11} & b_{12} & \cdots & b_{1p} \\ a_{21} & a_{22} & \cdots & a_{2n} & b_{21} & b_{22} & \cdots &b_{2p} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_{m1} & b_{m2} & \cdots & b_{mp} \end{bmatrix}$
当rank(A) = rank(B) = r < n,方程组有无穷多解
通解: $x = a_1x_1 + a_2 x_2 + \cdots + a_{n-r}x_{n-r} + x_0$

线性方程组的matlab求解
齐次方程组的通解
$\hat x = a_1 x_1 + a_2 x_2 + \cdots + a_{n-r} x_{n-r}$
求解A矩阵的化零矩阵, 使得 AZ = 0
Z = null(A)
特解: x_0 = pinv(A)*B

例4-43 方程求解
求解下列方程组的解
$\begin{bmatrix}1 & 4 & 0 & -1 & 0 & 7 & -9 \\ 2 & 8 & -1& 3 & 9 & -13 & 7 \\ 0 & 0 & 2 & -3& -4 & 12 & -8 \\ -1 & -4 & 2 & 4 & 8 & -31 & 37 \end{bmatrix} x = \begin{bmatrix}3 \\ 9 \\ 1 \\ 4 \end{bmatrix}$
判定矩阵方程的可解性

A = [1, 4, 0, -1, 0, 7, -9; 2, 8, -1, 3, 9, -13, 7; 0, 0, 2, -3, -4, 12, -8; -1, -4, 2, 4, 8, -31, 37];
B = [3;9;1;4];
C = [A B];
ra = rank(A);
rc = rank(C);
ra, rc

两种解法

% 方法1：通解与检验
Z = null(sym(A)),
x0 = sym(pinv(A)*B)
syms a1 a2 a3 a4;
x = Z*[a1;a2;a3;a4] + x0,
E = A*x -B
% 方法2：基本行变换方法
C = [A B];
D = rref(C)

$自变量x_2,x_5,x_6,x_7,记作 b_1,b_2,b_3,b_4$
方程的解
$\left \{ \begin{array}{l}x_1 = -4b_1 -2b_2 - b_3 + 3b_4 + 4 \\ x_3 = -b_2 + 3b_3 - 5b_2 + 2 \\ x_4 = -2b_2 + 6b_3 -6b_4 + 1 \end{array} \right.$

例4-44 解方法
试求解线性代数方程组
$\begin{bmatrix}4 & 7 & 1 & 4 \\ 3 & 7 & 4 & 6 \end{bmatrix} x = \begin{bmatrix} 3 \\ 4 \end{bmatrix}$

A = [4, 7, 1, 4; 3, 7, 4, 6];
B = [3; 4];
C = [A B];
rank(A), rank(C)
syms a1 a2 b1 b2;
x1 = null(sym(A))*[a1; a2] + sym(A \ B),
A*x1 - B,
a = rref(sym([A B]))
x2 = [a(:, 3:5)*[-b1; -b2; 1]; b1; b2],
A*x2 - B

III.方程组无解
如果rank(A) < rank(C),矛盾方程，无解
只能利用Moore-Penrose广义逆解出方程的最小二乘解
x = pinv(A)*B
它使误差的范围测度取最小值
$\min \limits_{x} \Vert Ax - B \Vert$

例4-45 矛盾方程
试求方程组
$\begin{bmatrix}1 & 2 & 3 & 4 \\ 2 & 2 & 1 & 1 \\ 2 & 4 & 6 & 8 \\ 4 & 4 &2 &2 \end{bmatrix} x = \begin{bmatrix}1 \\ 2 \\ 3 \\ 4 \end{bmatrix}$

A = [1, 2, 3, 4; 2, 2, 1, 1; 2, 4, 6, 8; 4, 4, 2, 2];
B = [1:4]';
C = [A B];
rank(A), rank(C)
% 最小二乘解
x = pinv(A)*B,
norm(A*x - B)

04.08 Lyapunov方程

离散Lyapunov方程
数学描述: $AXA^T - X + Q = 0$
函数调用格式 X = dlyap(A, Q)
Stein方程的特例: $B = A^T$
$(I_{n^2 \times n^2} - A \otimes A) x = q$

例4-45 解方程
求离散 Lyapunov 方程
$\begin{bmatrix}8 &1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2 \end{bmatrix} X \begin{bmatrix}8 &1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2 \end{bmatrix}^T - X + \begin{bmatrix} 16 & 4 & 1\\ 9 & 3 & 1 \\ 4 & 2 & 1 \end{bmatrix} = 0$

A = [8, 1, 6; 3, 5, 7; 4, 9, 2];
Q = [16, 4, 1; 9, 3, 1; 4, 2, 1];
X = dlyap(A, Q),
norm(A*X*A' - X + Q)

连续的Lyapunov方程
数学描述: $AX + XA^T = -C$
某些领域要求-C为对称正定的 n x n 矩阵
实际应用中 -C 还可以是任意矩阵
函数调用格式，该函数为控制系统工具想中的函数
X = lyap(A, C)

例4-46 Lyapunov方程
已知 $AX + XA^T = -C$
$A = \begin{bmatrix}1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 0 \end{bmatrix}, C = -\begin{bmatrix}10 & 5 & 4\\ 5 & 6 & 7 \\ 4 & 7 & 9 \end{bmatrix}$
求解相应的Lyapunov方程，并验证精确度

A = [1 2 3; 4 5 6 ; 7 8 0];
C = -[10 5 4; 5 6 7; 4 7 9];
X = lyap(A, C), 
norm(A*X+X*A' + C)

Lyapunov方程的解析解
方程的解 $(A \otimes I + I \otimes A) x = -c$
将Lyapunov方程的各个矩阵参数表示为
$X = \begin{bmatrix}x_1 & x_2 & \cdots & x_m \\ x_{m+1} & x_{m+2} & \cdots x_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ x_{(n-1)m+1} & x_{(n-1)m + 2} & \cdots & x_{nm} \end{bmatrix}$
$C = \begin{bmatrix}c_1 & c_2 & \cdots & c_m \\ c_{m+1} & c_{m+2} & \cdots & c_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ c_{(n-1)m+1} & c_{(m-1)m+2} & \cdots & c_{nm} \end{bmatrix}$

Kronecker乘积
$A \otimes B$, 为矩阵 A 和 B 的 Kronecker 乘积
$A \otimes B = \begin{bmatrix}a_{11}B & \cdots a_{1m}B \\ \vdots & \ddots & \vdots \\ a_{n1}B & \cdots & a_{nm}B \end{bmatrix}$
函数调用格式: C = kron(A, B)
一般情况下: $A \otimes B \neq B \otimes A$

例4-47 解析解
给定 $A = \begin{bmatrix}1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 0 \end{bmatrix}, C = \begin{bmatrix}10 & 5& 4 \\ 5 & 6 & 7 \\ 4 & 7 & 9 \end{bmatrix}$

试求解得到Lyapunov方程的解析解。

A0 = sym(kron(A,eye(3)) + kron(eye(3), A));
c = reshape(C', 9, 1);
x0 = -inv(A0) * c;
x = reshape(x0, 3, 3)'
norm(A*x + x*A' + C)

例4-48 复数方程
假定 $A = \begin{bmatrix}1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 0 \end{bmatrix}, C = -\begin{bmatrix}1+1j & 3 + 3j & 12 + 10j \\ 2 + 5j & 6 & 11+ 6j \\ 5 + 2j & 11 + j & 2 + 12j \end{bmatrix}$
c不为实对称正定矩阵的方程求解

A = [1 2 3; 4 5 6; 7 8 0];
C = -[1+1i, 3+3i, 12+10i; 2+5i, 6, 11+6i; 5+2i, 11+i, 2+12i];
A0 = sym(kron(A, eye(3)) + kron(eye(3), A));
c = reshape(C.', 9, 1);
x0 = -inv(A0)*c;
x = reshape(x0, 3, 3).',
norm(A*x+x*A.' + C)

Stein方程求解
数学描述: $AXB - X + Q = 0$
所有矩阵都为 n x n 方阵
令 x 为X矩阵的向量展开，q为Q矩阵的向量展开，Stein方程可以由下面的线性方程直接解出
$(I_{n^2 \times n^2} - B^T \otimes A) x = q$

例4-49 Stein 方程
求解Stein 方程
$\begin{bmatrix} -2 & 2 & 1 \\ -1 & 0 & -1 \\ 1 & -1 & 2 \end{bmatrix} X \begin{bmatrix} -2 & -1 & 2 \\ 1 & 3 & 0 \\ 3 & -2 & 2 \end{bmatrix} - X + \begin{bmatrix} 0 & 1 & 0 \\ -1 & 1 & 0 \\ 1 & -1 & -1 \end{bmatrix} = 0$

A = [-2, 2, 1; -1, 0, -1; 1, -1, 2];
B = [-2, -1, 2; 1, 3, 0; 3, -2, 2];
Q = [0, -1, 0; -1, 1, 0; 1, -1, -1];
x = inv(sym(eye(9)) - kron(B', A)) * Q(:);
X = reshape(x, 3, 3),
norm(A*X*B -X + Q)

04.09 Sylvester方程和Riccati方程

Sylvester方程求解
广义Lyapunov方程的数学描述
$AX + XB = -C$
函数调用格式 X = lyap(A, B, C)
使用Kronecker乘积得到解析解
$(A \otimes I_m + I_m \otimes B^T) x = c$

Sylvester方程的解析解
解析解:$(A \otimes I_m + I_m \otimes B^T) x = c$

function X = lyapsym(A, B, C)
  if nargin == 2
    C = B;
    B = A';
  end;
  [nr, nc] = size(C);
  A0 = kron(A, eye(nc)) + kron(eye(nr), B.');
  try
    C1 = C.';
    x0 = -inv(A0)*C1(:);
    X = reshape(x0, nc, nr).';
  catch
    error('singular matrix found.')
  end;
end

函数lyapsym()使用
连续Lyapunov方程
X = lyapsym(sym(A), C)
离散Lyapunov方程， 重新写成
$AX + X[-(A^T)^{-1}] = - Q(A^T)^{-1}$
X = lyapsym(sym(A), -inv(A.’), Q*inv(A.’))
Sylvester方程
X = lyapsym(sym(A), B, C)

例4-51 两种方法
解如下Sylvester方程
$\begin{bmatrix}8 & 1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2 \end{bmatrix} X + X \begin{bmatrix}16 & 4 & 1 \\ 9 & 3 & 1 \\ 4 & 2 & 1\end{bmatrix} = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 0 \end{bmatrix}$

A = [8, 1, 6; 3, 5, 7; 4, 9, 2];
B = [16, 4, 1; 9, 3, 1; 4, 2, 1];
C = -[1, 2, 3; 4, 5, 6; 7, 8, 0];
% 数值解法
X = lyap(A, B, C),
norm(A*X + X*B + C)
% 解析解法
x = lyapsym(sym(A), B, C),
norm(A*x + x*B + C)

例4-52 解析解
离散Lyapunov方程
$\begin{bmatrix}8 & 1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2 \end{bmatrix} X \begin{bmatrix}8 & 1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2 \end{bmatrix}^T - X + \begin{bmatrix}16 & 4 & 1\\ 9 & 3 & 1\\ 4 & 2 & 1 \end{bmatrix} = 0$

A = [8, 1, 6; 3, 5, 7; 4, 9, 2];
Q = [16, 4, 1; 9, 3, 1; 4, 2, 1];
x = lyapsym(sym(A), -inv(A.'), Q*inv(A.')),
norm(A*x*A.' - x + Q)

例4-53 长方形方程
求下面Sylvester方程
$A = \begin{bmatrix}8 & 1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2 \end{bmatrix}, B = \begin{bmatrix} 2 & 3 \\ 4 & 5 \end{bmatrix}, C = \begin{bmatrix}1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix}$

A = [8, 1, 6; 3, 5, 7; 4, 9, 2];
B = [2, 3; 4, 5];
C = [1, 2; 3, 4; 5, 6];
X = lyapsym(sym(A), B, C),
norm(A*X+X*B + C)

Riccati方程求解
Riccati方程数学描述: $A^TX + XA - XBX + C = 0$
二次型方程，没有解析解
函数调用格式: X = are(A, B, C)
注意：函数are()是基于数值运算的方法，ARE表示”algebraic Riccati equation.”

例4-54 解方程
$Riccati方程 A^TX + XA - XBX + C = 0$
$A = \begin{bmatrix}-2 & 1 & -3 \\ -1 & 0 & -2 \\ 0 & -1 & 2 \end{bmatrix}, B = \begin{bmatrix}2 & 2 & -2 \\ -1 & 5 & -2 \\ -1 & 1 & 2 \end{bmatrix}, C = \begin{bmatrix} 5 & -4 & 4 \\ 1 & 0 & 4 \\ 1 & -1 & 5 \end{bmatrix}$

A = [-2, 1, -3; -1, 0, -2; 0, -1, 2];
B = [2, 2, -2; -1, 5, -2; -1, 1, 2];
C = [5, -4, 4; 1, 0, 4; 1, -1, 5];
X = are(A, B, C),
norm(A'*X + X*A - X*B*X + C)

遗留问题
现有的函数are()只能求出方程的一个解
Riccati方程有其它解吗？
Riccati方程有多少解？有多少复数解？
如果有变形的Riccati方程，有办法求解吗？
$AX + XD - XBX + C = 0$
$AX + XD - XBX^T + C = 0$
第六章会介绍。

04.10 矩阵指数与三角函数

非线性运算与矩阵函数求值
矩阵的两类运算方法
面向矩阵元素的非线性预算、矩阵元素单独运算，相当于点运算
矩阵函数求值，这个矩阵运算

面向矩阵元素的非线性运算

函数名	运算意义
abs()	求模(绝对值)
asin(),acos(), atan()	反正弦，反余弦，反正切
sqrt()	求平方根
log(),log10(),exp()	自然对数、常用对数、指数函数
real(), imag(), conj()	求实部、虚部、共轭复数
sin(), cos(), tan()	正弦、余弦、正切
round(), floor(), ceil(), fix()	取整函数

例4-56 矩阵元素运算
$A = \begin{bmatrix}16 & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 5 & 1 \end{bmatrix}$
对A矩阵进行矩阵元素的指数和正弦运算

A = [16, 2, 3, 13; 5, 11, 10, 8; 9, 7, 6, 12; 4, 14, 15, 1];
exp(A),
sin(A)

矩阵函数求值
矩阵函数是对整个矩阵的函数运算，不再是对单独矩阵运算的运算

矩阵指数运算
Cleve Moler 文章总结了19种不同的数值算法，可以求解出矩阵指数 $e^A$
数值运算： E = expm(A)
符号运算：E = expm(sym(A))
矩阵的指数函数: F = expm(A*t)

例4-57 矩阵指数
$A = \begin{bmatrix}-2 & 1 & 0 & & \\ 0 & -2 & 1 & & \\ 0 & 0 & -2 & & \\ & & & -5 & 1 \\ & & & 0 & -5 \end{bmatrix}$

function A = diagm(varargin)
  A = [];
  for i = 1:length(varargin)
    A1 = varargin{i};
    [n, m] = size(A);
    [n1, m1] = size(A1);
    A(n+1:n+n1, m+1:m+m1) = A1;
  end
end

A1 = [-2, 1, 0; 0, -2, 1; 0, 0, -2];
A2 = [-5, 1; 0, -5];
A = diagm(A1, A2)
B = expm(A),
C = logm(B),
norm(C - A)
% 指数函数的解析解
syms t;
expm(A*t)

矩阵三角函数运算
matlab没有提供对矩阵进行三角函数的现成函数
求解任意非线性矩阵函数的函数调用格式
$A_1$ = funm(A, ‘fun_name’)
局限性：
采用的算法是矩阵的对角化方法，不能含有重特征根，否则不能得出正确结果
不能求出矩阵的三角函数 sin At

例4-59 矩阵的正弦
$A = \begin{bmatrix}-3 & -1 & -1 \\ 0 & -3 & -1 \\ 1 & 2 & 0 \end{bmatrix}$
求sin A

A = [-3, -1, -1; 0, -3, -1; 1, 2, 0];
B = funm(A, 'sin')

* 矩阵三角函数的解析解*
根据著名的Euler公式
$e^{ja} = \cos a + j \sin a$
$e^{-ja} = \cos a - j \sin a$
推导出
$\cos a = \dfrac{1}{2}(e^{ja} + e^{-ja})$
$\sin a = \dfrac{1}{j2}(e^{ja} - e^{-ja})$
此公式可以直接用于矩阵的形式

例4-68 矩阵三角函数的解析解
$A = \begin{bmatrix}-7 & 2 & 0 & -1 \\ 1 & -4 & 2 & 1 \\ 2 & -1 & -6 & -1 \\ -1 & -1 & 0 & 4 \end{bmatrix}$
求sin At 和 cos At

% 计算很慢
A = [-7, 2, 0, -1; 1, -4, 2, 1; 2, -1, -6, -1; -1, -1, 0, 4];
syms t;
A1 = (expm(A*1j*t) - expm(-A*1j*t))/(2*1j);
A2 = (expm(A*1j*t)+expm(-A*1j*t))/2;
A1 = simplify(A1),
A2 = simplify(A2)
simplify(funm(A*t, 'sin')),
simplify(funm(A*t, 'cos'))

04.11 矩阵任意函数计算

例 4-69 幂零矩阵的性质
幂零矩阵(nilpotent matrix)
$H = \begin{bmatrix}0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{bmatrix}$
幂零矩阵性质

H = diag([1 1 1], 1);
for i = 2: 4
  H^i
end

幂零矩阵的高次方是零矩阵

一般矩阵函数的运算
将 $m_i \times m_i$ Jordan 块记成 $J_i = \lambda_i I_{m_i} + H_{m_i \times m_i}$
其中， $\lambda_i$ 为 Jordan 矩阵的 $m_i$ 重特征值
$H_{m_i}$ 为幂零矩阵，即下式成立
$k \geq m_i, H_{m_i}^k \equiv 0$
矩阵函数
$\psi(J_i) = \psi(\lambda_i) I_{m_i} + \cdots + \dfrac{\psi^{(m_i - 1)}(\lambda_i)}{(m_i-1)!}H_{m_i}^{m_i-1}$
将任意已知矩阵A进行Jordan分解
$A = V \begin{bmatrix}J_1 & & & \\ & J_2 & & \\ & & \ddots & \\ & & & J_m \end{bmatrix} V^{-1}$
则对该矩阵的任意函数运算符如下
$\psi(A) = V \begin{bmatrix}\psi(J_1) & & & \\ & \psi(J_2) & & \\ & & \ddots & \\ & & & \psi(J_m) \end{bmatrix}V^{-1}$

% 求矩阵函数的解析解
function F = funmsym(A, fun, x)
  [V,T] = jordan(A);
  vec = diag(T);
  v1 = [0,diag(T,1)', 0]; 
  v2 = find(v1 == 0);
  lam = vec(v2(1: end-1));
  m = length(lam); 
  for i = 1:m
    k = v2(i):v2(i+1)-1;
    J1 = T(k, k);
    F(k, k) = funJ(J1, fun, x);
  end
  F = V*F*inv(V);
function fJ = funJ(J, fun, x)
  lam = J(1,1);
  f1 = fun;
  fJ = subs(fun, x, lam)*eye(size(J));
  H = diag(diag(J,1),1);
  H1 = H;
  for i = 2:length(J)
    f1 = diff(f1,x);
    a1 = subs(f1,x,lam);
    fJ = fJ+a1*H1;
    H1 = H1*H/i;
  end

例 4-63 复合矩阵函数
$A = \begin{bmatrix}-7 & 2 & 0 & -1 \\ 1 & -4 & 2 & 1 \\ 2 & -1 & -6 & -1 \\ -1 & -1 & 0 & -4 \end{bmatrix}$
求下列矩阵函数 $\psi(A) = e^{A \cos At}$

A = [-7, 2, 0, -1; 1, -4, 2, 1; 2, -1, -6, -1; -1, -1, 0, -4];
syms x t;
A1 = funmsym(A, exp(x*cos(x*t)), x)
collect(A1(1, 1), exp(-6*cos(6*t)))
subs(A1, t, 1)

矩阵乘方计算 – $A^k$
如何结算$A^k$?
仿照矩阵任意函数求解，利用Jordan分解与幂零矩阵
如果 $A = VJV^{-1}, 则 A^k = VJ^kV^{-1}$
如果 $J = \lambda I + H_m$
$J^k = \lambda^k I + k \lambda^{k-1} H_m + \dfrac{k(k-1)}{2!} \lambda^{k-2} H_m^2 + \cdots$

例4-71 矩阵乘方计算
$A = \begin{bmatrix}-3 & -1 & -1 \\ 0 & -3 & -1 \\ 1 & 2 & 0 \end{bmatrix}$

% 乘方计算
A = sym([-3, -1, -1; 0, -3, -1; 1, 2, 0]);
syms k;
[V J] = jordan(sym(A))
A0 = 2*eye(3);
H = J - A0;
J1 = A0^k + k*A0^(k-1) * H + k*(k-1)/2*A0^(k-2)*H^2;
F = simplify(V*J1*inv(V))

矩阵乘方$A^k$自定义函数

function F = mpowersym(A, k)
  A = sym(A);
  [V,T] = jordan(A);
  vec = diag(T);
  v1 = [0, diag(T,1)', 0];
  v2 = find(v1 == 0);
  lam = vec(v2(1:end-1));
  m = length(lam);
  for i = 1:m,
    k0 = v2(i):v2(i+1)-1;
    J1 = T(k0,k0);
    F(k0,k0) = powJ(J1,k); 
  end
  F = simplify(V*F*inv(V));

function fJ = powJ(J,k)
  lam = J(1,1);
  I = eye(size(J));
  H = J-lam*I;
  fJ = lam^k*I;
  H1 = k*H; 
  for i = 2:length(J)
    fJ = fJ + lam^(k+1-i)*I*H1;
    H1 = H1*H*(k+1-i)/i;
  end

例4-72 矩阵乘方计算
$A = \begin{bmatrix}-7 & 2 & 0 & -1 \\ 1 & -4 & 2 & 1 \\ 2 & -1 & -6 & -1 \\ -1 & -1 & 0 & -4 \end{bmatrix}$

A = [-7, 2, 0, -1; 1, -4, 2, 1; 2, -1, -6, -1; -1, -1, 0, -4];
syms k;
A = sym(A);
F = mpowersym(A, k)

矩阵乘方检验
检验A的平方根与$A^{12345}$

simplify(A^12345 - subs(F, k, 12345))
syms x;
A3 = funmsym(sym(A), sqrt(x), x),
simplify(A3 - subs(F, k, 1/2))