现代科学运算-matlab语言与应用
东北大学 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(μ,σ2):B=μ+σ∗randn(n,m)
N
(
μ
,
σ
2
)
:
B
=
μ
+
σ
∗
r
a
n
d
n
(
n
,
m
)
对角元素矩阵
对角矩阵的数学描述
diag(a1,a2,⋯,an)=⎡⎣⎢⎢⎢⎢⎢a1a2⋱an⎤⎦⎥⎥⎥⎥⎥
d
i
a
g
(
a
1
,
a
2
,
⋯
,
a
n
)
=
[
a
1
a
2
⋱
a
n
]
其中,所有的非对角元素都为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=⎡⎣⎢⎢⎢⎢⎢c1c2⋮cnc2c3⋮cn+1⋯⋯⋱⋯cmcm+1⋮cn+m−1⎤⎦⎥⎥⎥⎥⎥
H
=
[
c
1
c
2
⋯
c
m
c
2
c
3
⋯
c
m
+
1
⋮
⋮
⋱
⋮
c
n
c
n
+
1
⋯
c
n
+
m
−
1
]
反对角线上元素相同
函数调用:
H = hankel(C)
H = hankel(C, R)
C(end) = R(1)
Vandermonde矩阵
Vandermonde矩阵的数学描述
V=⎡⎣⎢⎢⎢⎢⎢cn−11cn−12⋮cn−1ncn−22cn−22⋮cn−2n⋯⋯⋱⋯c1c2⋮cn11⋮1⎤⎦⎥⎥⎥⎥⎥
V
=
[
c
1
n
−
1
c
2
n
−
2
⋯
c
1
1
c
2
n
−
1
c
2
n
−
2
⋯
c
2
1
⋮
⋮
⋱
⋮
⋮
c
n
n
−
1
c
n
n
−
2
⋯
c
n
1
]
其中,
vi,j=cn−ji,i,j=1,2,⋯,n
v
i
,
j
=
c
i
n
−
j
,
i
,
j
=
1
,
2
,
⋯
,
n
生成Vandermonde矩阵 V = vander(c)
相伴矩阵
一个首一化的多项式
P(s)=sn+a1sn−1+a2sn−2+⋯+an−1s+an
P
(
s
)
=
s
n
+
a
1
s
n
−
1
+
a
2
s
n
−
2
+
⋯
+
a
n
−
1
s
+
a
n
其相伴矩阵
Ac=⎡⎣⎢⎢⎢⎢⎢⎢⎢−a110⋮0−a201⋮0⋯⋯⋯⋱⋯−an−100⋮1−an00⋮0⎤⎦⎥⎥⎥⎥⎥⎥⎥
A
c
=
[
−
a
1
−
a
2
⋯
−
a
n
−
1
−
a
n
1
0
⋯
0
0
0
1
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
1
0
]
matlab函数: B = compan(p)
例4-5 由多项式构造相伴矩阵
考虑一个多项式
P(s)=2s4+4s2+5s+6
P
(
s
)
=
2
s
4
+
4
s
2
+
5
s
+
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=⎡⎣⎢⎢⎢a11a21a31a41a12a22a32a42a13a23a33a43a14a24a34a44⎤⎦⎥⎥⎥
A
=
[
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
]
B=⎡⎣⎢⎢⎢a11a21a31a41a12a22a32a42⎤⎦⎥⎥⎥
B
=
[
a
11
a
12
a
21
a
22
a
31
a
32
a
41
a
42
]
C=⎡⎣⎢⎢⎢f11f21f31f41f12f22f32f42f13f23f33f43f14f24f34f44⎤⎦⎥⎥⎥
C
=
[
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
]
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={aij}
A
=
{
a
i
j
}
的行列式定义为
D=|A|=det(A)=∑(−1)ka1k1a2k2⋯ankn
D
=
|
A
|
=
d
e
t
(
A
)
=
∑
(
−
1
)
k
a
1
k
1
a
2
k
2
⋯
a
n
k
n
函数调用:
d = det(A), d = det(sym(A))
注意: 该函数可用于数值和符号运算,方法和A的类型一致
自动识别矩阵的形式,选择求解方法
例4-8 矩阵求行列式
A=⎡⎣⎢⎢⎢16594211714310615138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
15
1
]
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 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={aij},i,j=1,2,⋯,n
A
=
{
a
i
j
}
,
i
,
j
=
1
,
2
,
⋯
,
n
矩阵的迹定义为对角线元素之和:tr(A) =
∑i=1naij
∑
i
=
1
n
a
i
j
函数调用格式:t = trace(A)
迹也等于特征值的和:sum(diag(A))
矩阵的秩
矩阵A的秩:线性不相关的最大行(列)数
rank(A)
=rc=rr
=
r
c
=
r
r
数值方法或符号方法: r = rank(A)
给定精度
ε
ε
下求数值秩: r = rank(A,
ε
ε
)
向量的范数
范数式测度 – 用一个值描述向量、矩阵的大小
函数ρ(x)为x向量的范数的条件:
函
数
ρ
(
x
)
为
x
向
量
的
范
数
的
条
件
:
ρ(x)≥0且ρ(x)=0的充要条件是x=0
ρ
(
x
)
≥
0
且
ρ
(
x
)
=
0
的
充
要
条
件
是
x
=
0
ρ(ax)=|a|ρ(x),a为任意标量
ρ
(
a
x
)
=
|
a
|
ρ
(
x
)
,
a
为
任
意
标
量
对向量x和y,有ρ(x+y)≤ρ(x)+ρ(y)
对
向
量
x
和
y
,
有
ρ
(
x
+
y
)
≤
ρ
(
x
)
+
ρ
(
y
)
向量的p-范数
∥x∥∞=max1≤i≤n|xi|
‖
x
‖
∞
=
max
1
≤
i
≤
n
|
x
i
|
∥x∥p=(∑i=1n|xi|p)1p,p=1,2,⋯,
‖
x
‖
p
=
(
∑
i
=
1
n
|
x
i
|
p
)
1
p
,
p
=
1
,
2
,
⋯
,
矩阵的范数
对于任意非零向量x, 矩阵A的范数是
∥A∥=supx≠0∥Ax∥∥x∥
‖
A
‖
=
sup
x
≠
0
‖
A
x
‖
‖
x
‖
其它常用范数
∥A∥1=max1≤j≤n∑i=1n|aij|,∥A∥2=smax(ATA)−−−−−−−−−√
‖
A
‖
1
=
max
1
≤
j
≤
n
∑
i
=
1
n
|
a
i
j
|
,
‖
A
‖
2
=
s
m
a
x
(
A
T
A
)
∥A∥∞=max1≤i≤n∑i=1n|aij|,∥A∥F=trace(ATA)
‖
A
‖
∞
=
max
1
≤
i
≤
n
∑
i
=
1
n
|
a
i
j
|
,
‖
A
‖
F
=
t
r
a
c
e
(
A
T
A
)
矩阵范数计算
s(X)为x矩阵的特征值,
smax(ATA)
s
m
a
x
(
A
T
A
)
A矩阵的范数为矩阵AA^T其最大特征值
函数调用格式:
∥A∥2
‖
A
‖
2
N = norm(A)
选项可以是1,2, inf, ‘fro’: N = norm(A, opts)
注意: 早期版本norm()函数仅用于数值矩阵
矩阵特征多项式
由矩阵A的如下多项式:
c(s)=det(sI−A)=sn+c1sn−1+⋯+cn−1s+cn
c
(
s
)
=
d
e
t
(
s
I
−
A
)
=
s
n
+
c
1
s
n
−
1
+
⋯
+
c
n
−
1
s
+
c
n
多项式c(s)是矩阵A的特征多项式
matlab下多项式由降幂排列的系数向量表示
c = poly(A)
c = charpoly(sym(A))
c = charpoly(sym(A), x)
矩阵多项式求解
矩阵多项式的数学表示:
B=a1An+a2An−1+⋯+anA+an+1I
B
=
a
1
A
n
+
a
2
A
n
−
1
+
⋯
+
a
n
A
+
a
n
+
1
I
函数调用格式,A不能为符号矩阵
B = polyvalm(a, A)
其中,a = [
a1,a2,⋯,an,an+1
a
1
,
a
2
,
⋯
,
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 =
a1
a
1
x .^n +
a2
a
2
x.^(n-1) + … +
an+1
a
n
+
1
调用格式:C = polyval(a, x)
给出多项式p(符号运算工具箱):C = subs(p, s, x)
Cayley-Hamilton定理
若矩阵A的特征多项式为
f(x)=det(sI−A)=a1sn+a2sn−1+⋯+ans+an+1
f
(
x
)
=
d
e
t
(
s
I
−
A
)
=
a
1
s
n
+
a
2
s
n
−
1
+
⋯
+
a
n
s
+
a
n
+
1
则
f(A)=0
f
(
A
)
=
0
, 即
a1An+a2An−1+⋯+anA+an+1I=0
a
1
A
n
+
a
2
A
n
−
1
+
⋯
+
a
n
A
+
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
=
A
−
1
函数表示: C = inv(A)
基本行变换(行阶梯型矩阵):
H1=rref(H)
H
1
=
r
r
e
f
(
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=⎡⎣⎢⎢⎢16594211714310615138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
15
1
]
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=⎡⎣⎢⎢⎢a1a2a3a4a2a3a40a3a400a4000⎤⎦⎥⎥⎥
A
=
[
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
]
a = sym('a%d', [1, 4]);
H = hankel(a);
inv(H)
矩阵的广义逆
适用于奇异或矩形矩阵
如果存在ANA = A,则称 N为A矩阵的广义逆矩阵,记作:
N=A−
N
=
A
−
不唯一
定义范数最小化指标
minN∥AN−I∥
min
N
‖
A
N
−
I
‖
Moore-Penrose广义逆
矩阵M为矩阵A的Moore-Penrose广义逆矩阵的条件:
(i) AMA = A
(ii) MAM = M
(iii) AM 和 MA 均为 Hermite 对称矩阵
记为
M=A+
M
=
A
+
, 又称 伪逆
Moore-Penrose 广义逆矩阵是唯一的
M = pinv(A), M = pinv(A,
ϵ
ϵ
)
例4-24 奇异矩阵的伪逆
使用pinv()函数计算矩阵A的伪逆
A=⎡⎣⎢⎢⎢16594211714310615138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
15
1
]
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=⎡⎣⎢63−310−241−5248124⎤⎦⎥
A
=
[
6
1
4
2
1
3
0
1
4
2
−
3
−
2
−
5
8
4
]
% 求秩并求伪逆
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=λx
A
x
=
λ
x
非零向量x是特征向量,数值
λ
λ
是特征值
调用格式: d = eig(A) 或 [V, D] = eig(A)
例4-26 矩阵的特征值
计算下述矩阵的特征值与特征向量
A=⎡⎣⎢⎢⎢16594211714310615138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
15
1
]
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−1AV=D
V
−
1
A
V
=
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,c1,c2,c3]
c
=
[
1
,
c
1
,
c
2
,
c
3
]
可以构造出相伴矩阵
H=⎡⎣⎢−c110−c201−c300⎤⎦⎥
H
=
[
−
c
1
−
c
2
−
c
3
1
0
0
0
1
0
]
试求其特征值与特征向量矩阵
syms c1 c2 c3 real;
H = compan([1 c1 c2 c3]);
[V, D] = eig(H)
矩阵的广义特征向量问题
数学描述:
Ax=λBx
A
x
=
λ
B
x
非零向量x是特征向量, 数值
λ
λ
是特征值,
B 是正定对称矩阵
调用格式 d = eig(A, B)
或 [V, D] = eig(A, B)
例4-27 广义特征值
已知A, B,求广义特征值与特征向量矩阵
A=⎡⎣⎢⎢⎢5765710876810957910⎤⎦⎥⎥⎥,B=⎡⎣⎢⎢⎢25−356−1−4−2−121−3−23108⎤⎦⎥⎥⎥
A
=
[
5
7
6
5
7
10
8
7
6
8
10
9
5
7
9
10
]
,
B
=
[
2
6
−
1
−
2
5
−
1
2
3
−
3
−
4
1
10
5
−
2
−
3
8
]
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−1AB
X
=
B
−
1
A
B
正交矩阵的数学描述
QHQ=I,QQH=I
Q
H
Q
=
I
,
Q
Q
H
=
I
函数调用格式: Q = orth(A)
矩阵的三角分解
一般矩阵分解成下三角与上三角矩阵的积
数学描述: A = LU
其中
L=⎡⎣⎢⎢⎢⎢⎢1l21⋮ln11⋮ln2⋱⋯1⎤⎦⎥⎥⎥⎥⎥,U=⎡⎣⎢⎢⎢⎢⎢u11u12u22⋯⋯⋱u1nu2n⋮unn⎤⎦⎥⎥⎥⎥⎥
L
=
[
1
l
21
1
⋮
⋮
⋱
l
n
1
l
n
2
⋯
1
]
,
U
=
[
u
11
u
12
⋯
u
1
n
u
22
⋯
u
2
n
⋱
⋮
u
n
n
]
三角分解的递推算法
递推初值
u1i=a1i,i=1,2,⋯,n
u
1
i
=
a
1
i
,
i
=
1
,
2
,
⋯
,
n
递推计算公式
lij=aij−∑k=1j−1likukjujj,(j<i)
l
i
j
=
a
i
j
−
∑
k
=
1
j
−
1
l
i
k
u
k
j
u
j
j
,
(
j
<
i
)
uij=aij−∑k=1i−1likukj,(j≥i)
u
i
j
=
a
i
j
−
∑
k
=
1
i
−
1
l
i
k
u
k
j
,
(
j
≥
i
)
三角分解的matlab求解
lu分解:
数学描述:A = LU
matlab求解: [L, U] = lu(A)
P为置换矩阵
数学表示:
A=P−1LU
A
=
P
−
1
L
U
matlab求解: [L, U, P] = lu(A)
例4-30 矩阵的三角分解
给定
A=⎡⎣⎢⎢⎢16594211714310615138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
15
1
]
两种方法调用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=LLT
A
=
L
L
T
下三角矩阵:
L=⎡⎣⎢⎢⎢⎢⎢l11l21⋮ln1l22⋮ln2⋱⋯lnn⎤⎦⎥⎥⎥⎥⎥
L
=
[
l
11
l
21
l
22
⋮
⋮
⋱
l
n
1
l
n
2
⋯
l
n
n
]
其中 A 是对称矩阵
Cholesky分解的算法
对称矩阵的Cholesky分解算法
lii=aii−∑k=1i−1l2ik−−−−−−−−−√
l
i
i
=
a
i
i
−
∑
k
=
1
i
−
1
l
i
k
2
lji=1ljj(aij−∑k=1j−1likljk),j<i
l
j
i
=
1
l
j
j
(
a
i
j
−
∑
k
=
1
j
−
1
l
i
k
l
j
k
)
,
j
<
i
调用格式: D = chol(A)
例4-26 Cholesky分解
对如下A矩阵进行Cholesky分解
A=⎡⎣⎢⎢⎢9342360740602709⎤⎦⎥⎥⎥
A
=
[
9
3
4
2
3
6
0
7
4
0
6
0
2
7
0
9
]
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,⋯,An−1x]
T
=
[
x
,
A
x
,
⋯
,
A
n
−
1
x
]
矩阵A可以转换成一个与类似相伴矩阵的形式
这样的x向量有无穷多个
例4-33 寻找变换矩阵
A=⎡⎣⎢⎢⎢5765710876810957910⎤⎦⎥⎥⎥
A
=
[
5
7
6
5
7
10
8
7
6
8
10
9
5
7
9
10
]
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=⎡⎣⎢⎢⎢31−1022−212−2032−2−25⎤⎦⎥⎥⎥
A
=
[
3
2
2
2
1
2
−
2
−
2
−
1
−
2
0
−
2
0
1
3
5
]
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=⎡⎣⎢⎢⎢−71750−67−65894−121−811178−173−46504−58⎤⎦⎥⎥⎥
A
=
[
−
71
−
65
−
81
−
46
75
89
117
50
0
4
8
4
−
67
−
121
−
173
−
58
]
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))
矩阵的奇异值分解
数学描述:
ATA≥0,AAT≥0
A
T
A
≥
0
,
A
A
T
≥
0
理论上,
rank(ATA)=rank(AAT)=rank(A)
r
a
n
k
(
A
T
A
)
=
r
a
n
k
(
A
A
T
)
=
r
a
n
k
(
A
)
奇异值定义:
σi(A)=λi(ATA)−−−−−−−√
σ
i
(
A
)
=
λ
i
(
A
T
A
)
其中,
σi
σ
i
是非负特征值
例4-40 奇异值实例
A=⎡⎣⎢1μ010μ⎤⎦⎥μ=10−15
A
=
[
1
1
μ
0
0
μ
]
μ
=
10
−
15
求矩阵的秩
ATA=[1+μ2111+μ2]
A
T
A
=
[
1
+
μ
2
1
1
1
+
μ
2
]
A = [1, 1; 5*eps, 0; 0, 5*eps];
rank(A)
奇异值分解
奇异值分解 Singular value decomposition(SVD)
A=LA1M
A
=
L
A
1
M
A – n x m 矩阵
L 和 M 为正交矩阵
分解出对角矩阵
A1=diag(σ1,⋯,σn)
A
1
=
d
i
a
g
(
σ
1
,
⋯
,
σ
n
)
对角元素满足:
σ1≥σ2≥⋯≥σn≥0
σ
1
≥
σ
2
≥
⋯
≥
σ
n
≥
0
求解: S = svd(A)
[L,
A1
A
1
, M] = svd(A)
矩阵的条件数
矩阵的条件数
cond(A)=σmaxσmin=σ¯(A)σ––(A)
c
o
n
d
(
A
)
=
σ
max
σ
min
=
σ
¯
(
A
)
σ
_
(
A
)
matlab求解: c = cond(A)
条件数过大,说明矩阵接近奇异值
坏条件矩阵
处理时应该慎重(建议用符号运算)
例4-40 奇异值分解
对下列矩阵进行奇异值分解
A=⎡⎣⎢⎢⎢16594211714310615138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
15
1
]
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=[12345678]
A
=
[
1
3
5
7
2
4
6
8
]
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=⎡⎣⎢⎢⎢⎢⎢a11a21⋮am1a12a22⋮am2⋯⋯⋱⋯a1na2n⋮amn⎤⎦⎥⎥⎥⎥⎥,B=⎡⎣⎢⎢⎢⎢⎢b11b21⋮bm1b12b22⋮bm2⋯⋯⋱⋯b1pb2p⋮bmp⎤⎦⎥⎥⎥⎥⎥
A
=
[
a
11
a
12
⋯
a
1
n
a
21
a
22
⋯
a
2
n
⋮
⋮
⋱
⋮
a
m
1
a
m
2
⋯
a
m
n
]
,
B
=
[
b
11
b
12
⋯
b
1
p
b
21
b
22
⋯
b
2
p
⋮
⋮
⋱
⋮
b
m
1
b
m
2
⋯
b
m
p
]
线性代数课程:解的三种情况
I.唯一解
唯一解存在条件
当 m = n,且 rank(A) = n, 则唯一解存在
X=A−1B
X
=
A
−
1
B
函数调用格式:
X = inv(A)*B
X = A \ B
注意: 推荐使用符号运算方法
例4-42 解方程
求解线性方程组
⎡⎣⎢⎢⎢1414233132234142⎤⎦⎥⎥⎥X=⎡⎣⎢⎢⎢54321234⎤⎦⎥⎥⎥
[
1
2
3
4
4
3
2
1
1
3
2
4
4
1
3
2
]
X
=
[
5
1
4
2
3
3
2
4
]
求解与检验
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=⎡⎣⎢⎢⎢⎢⎢a11a21⋮am1a12a22⋮am2⋯⋯⋱⋯a1na2n⋮amnb11b21⋮bm1b12b22⋮bm2⋯⋯⋱⋯b1pb2p⋮bmp⎤⎦⎥⎥⎥⎥⎥
C
=
[
a
11
a
12
⋯
a
1
n
b
11
b
12
⋯
b
1
p
a
21
a
22
⋯
a
2
n
b
21
b
22
⋯
b
2
p
⋮
⋮
⋱
⋮
⋮
⋮
⋱
⋮
a
m
1
a
m
2
⋯
a
m
n
b
m
1
b
m
2
⋯
b
m
p
]
当rank(A) = rank(B) = r < n,方程组有无穷多解
通解:
x=a1x1+a2x2+⋯+an−rxn−r+x0
x
=
a
1
x
1
+
a
2
x
2
+
⋯
+
a
n
−
r
x
n
−
r
+
x
0
线性方程组的matlab求解
齐次方程组的通解
x^=a1x1+a2x2+⋯+an−rxn−r
x
^
=
a
1
x
1
+
a
2
x
2
+
⋯
+
a
n
−
r
x
n
−
r
求解A矩阵的化零矩阵, 使得 AZ = 0
Z = null(A)
特解: x_0 = pinv(A)*B
例4-43 方程求解
求解下列方程组的解
⎡⎣⎢⎢⎢120−1480−40−122−13−3409−487−1312−31−97−837⎤⎦⎥⎥⎥x=⎡⎣⎢⎢⎢3914⎤⎦⎥⎥⎥
[
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
]
x
=
[
3
9
1
4
]
判定矩阵方程的可解性
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)
自变量x2,x5,x6,x7,记作b1,b2,b3,b4
自
变
量
x
2
,
x
5
,
x
6
,
x
7
,
记
作
b
1
,
b
2
,
b
3
,
b
4
方程的解
⎧⎩⎨⎪⎪x1=−4b1−2b2−b3+3b4+4x3=−b2+3b3−5b2+2x4=−2b2+6b3−6b4+1
{
x
1
=
−
4
b
1
−
2
b
2
−
b
3
+
3
b
4
+
4
x
3
=
−
b
2
+
3
b
3
−
5
b
2
+
2
x
4
=
−
2
b
2
+
6
b
3
−
6
b
4
+
1
例4-44 解方法
试求解线性代数方程组
[43771446]x=[34]
[
4
7
1
4
3
7
4
6
]
x
=
[
3
4
]
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
它使误差的范围测度取最小值
minx∥Ax−B∥
min
x
‖
A
x
−
B
‖
例4-45 矛盾方程
试求方程组
⎡⎣⎢⎢⎢1224224431624182⎤⎦⎥⎥⎥x=⎡⎣⎢⎢⎢1234⎤⎦⎥⎥⎥
[
1
2
3
4
2
2
1
1
2
4
6
8
4
4
2
2
]
x
=
[
1
2
3
4
]
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方程
数学描述:
AXAT−X+Q=0
A
X
A
T
−
X
+
Q
=
0
函数调用格式 X = dlyap(A, Q)
Stein方程的特例:
B=AT
B
=
A
T
(In2×n2−A⊗A)x=q
(
I
n
2
×
n
2
−
A
⊗
A
)
x
=
q
例4-45 解方程
求离散 Lyapunov 方程
⎡⎣⎢834159672⎤⎦⎥X⎡⎣⎢834159672⎤⎦⎥T−X+⎡⎣⎢1694432111⎤⎦⎥=0
[
8
1
6
3
5
7
4
9
2
]
X
[
8
1
6
3
5
7
4
9
2
]
T
−
X
+
[
16
4
1
9
3
1
4
2
1
]
=
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+XAT=−C
A
X
+
X
A
T
=
−
C
某些领域要求-C为对称正定的 n x n 矩阵
实际应用中 -C 还可以是任意矩阵
函数调用格式,该函数为控制系统工具想中的函数
X = lyap(A, C)
例4-46 Lyapunov方程
已知
AX+XAT=−C
A
X
+
X
A
T
=
−
C
A=⎡⎣⎢147258360⎤⎦⎥,C=−⎡⎣⎢1054567479⎤⎦⎥
A
=
[
1
2
3
4
5
6
7
8
0
]
,
C
=
−
[
10
5
4
5
6
7
4
7
9
]
求解相应的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⊗I+I⊗A)x=−c
(
A
⊗
I
+
I
⊗
A
)
x
=
−
c
将Lyapunov方程的各个矩阵参数表示为
X=⎡⎣⎢⎢⎢⎢⎢x1xm+1⋮x(n−1)m+1x2xm+2⋮x(n−1)m+2⋯⋯x2m⋱⋯xm⋮xnm⎤⎦⎥⎥⎥⎥⎥
X
=
[
x
1
x
2
⋯
x
m
x
m
+
1
x
m
+
2
⋯
x
2
m
⋮
⋮
⋱
⋮
x
(
n
−
1
)
m
+
1
x
(
n
−
1
)
m
+
2
⋯
x
n
m
]
C=⎡⎣⎢⎢⎢⎢⎢c1cm+1⋮c(n−1)m+1c2cm+2⋮c(m−1)m+2⋯⋯⋱⋯cmc2m⋮cnm⎤⎦⎥⎥⎥⎥⎥
C
=
[
c
1
c
2
⋯
c
m
c
m
+
1
c
m
+
2
⋯
c
2
m
⋮
⋮
⋱
⋮
c
(
n
−
1
)
m
+
1
c
(
m
−
1
)
m
+
2
⋯
c
n
m
]
Kronecker乘积
A⊗B
A
⊗
B
, 为矩阵 A 和 B 的 Kronecker 乘积
A⊗B=⎡⎣⎢⎢a11B⋮an1B⋯a1mB⋱⋯⋮anmB⎤⎦⎥⎥
A
⊗
B
=
[
a
11
B
⋯
a
1
m
B
⋮
⋱
⋮
a
n
1
B
⋯
a
n
m
B
]
函数调用格式: C = kron(A, B)
一般情况下:
A⊗B≠B⊗A
A
⊗
B
≠
B
⊗
A
例4-47 解析解
给定
A=⎡⎣⎢147258360⎤⎦⎥,C=⎡⎣⎢1054567479⎤⎦⎥
A
=
[
1
2
3
4
5
6
7
8
0
]
,
C
=
[
10
5
4
5
6
7
4
7
9
]
试求解得到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=⎡⎣⎢147258360⎤⎦⎥,C=−⎡⎣⎢1+1j2+5j5+2j3+3j611+j12+10j11+6j2+12j⎤⎦⎥
A
=
[
1
2
3
4
5
6
7
8
0
]
,
C
=
−
[
1
+
1
j
3
+
3
j
12
+
10
j
2
+
5
j
6
11
+
6
j
5
+
2
j
11
+
j
2
+
12
j
]
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
A
X
B
−
X
+
Q
=
0
所有矩阵都为 n x n 方阵
令 x 为X矩阵的向量展开,q为Q矩阵的向量展开,Stein方程可以由下面的线性方程直接解出
(In2×n2−BT⊗A)x=q
(
I
n
2
×
n
2
−
B
T
⊗
A
)
x
=
q
例4-49 Stein 方程
求解Stein 方程
⎡⎣⎢−2−1120−11−12⎤⎦⎥X⎡⎣⎢−213−13−2202⎤⎦⎥−X+⎡⎣⎢0−1111−100−1⎤⎦⎥=0
[
−
2
2
1
−
1
0
−
1
1
−
1
2
]
X
[
−
2
−
1
2
1
3
0
3
−
2
2
]
−
X
+
[
0
1
0
−
1
1
0
1
−
1
−
1
]
=
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
A
X
+
X
B
=
−
C
函数调用格式 X = lyap(A, B, C)
使用Kronecker乘积得到解析解
(A⊗Im+Im⊗BT)x=c
(
A
⊗
I
m
+
I
m
⊗
B
T
)
x
=
c
Sylvester方程的解析解
解析解:
(A⊗Im+Im⊗BT)x=c
(
A
⊗
I
m
+
I
m
⊗
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[−(AT)−1]=−Q(AT)−1
A
X
+
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方程
⎡⎣⎢834159672⎤⎦⎥X+X⎡⎣⎢1694432111⎤⎦⎥=⎡⎣⎢147258360⎤⎦⎥
[
8
1
6
3
5
7
4
9
2
]
X
+
X
[
16
4
1
9
3
1
4
2
1
]
=
[
1
2
3
4
5
6
7
8
0
]
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方程
⎡⎣⎢834159672⎤⎦⎥X⎡⎣⎢834159672⎤⎦⎥T−X+⎡⎣⎢1694432111⎤⎦⎥=0
[
8
1
6
3
5
7
4
9
2
]
X
[
8
1
6
3
5
7
4
9
2
]
T
−
X
+
[
16
4
1
9
3
1
4
2
1
]
=
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=⎡⎣⎢834159672⎤⎦⎥,B=[2435],C=⎡⎣⎢135246⎤⎦⎥
A
=
[
8
1
6
3
5
7
4
9
2
]
,
B
=
[
2
3
4
5
]
,
C
=
[
1
2
3
4
5
6
]
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方程数学描述:
ATX+XA−XBX+C=0
A
T
X
+
X
A
−
X
B
X
+
C
=
0
二次型方程,没有解析解
函数调用格式: X = are(A, B, C)
注意:函数are()是基于数值运算的方法,ARE表示”algebraic Riccati equation.”
例4-54 解方程
Riccati方程ATX+XA−XBX+C=0
R
i
c
c
a
t
i
方
程
A
T
X
+
X
A
−
X
B
X
+
C
=
0
A=⎡⎣⎢−2−1010−1−3−22⎤⎦⎥,B=⎡⎣⎢2−1−1251−2−22⎤⎦⎥,C=⎡⎣⎢511−40−1445⎤⎦⎥
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
]
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
A
X
+
X
D
−
X
B
X
+
C
=
0
AX+XD−XBXT+C=0
A
X
+
X
D
−
X
B
X
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=⎡⎣⎢⎢⎢1659421171431065138121⎤⎦⎥⎥⎥
A
=
[
16
2
3
13
5
11
10
8
9
7
6
12
4
14
5
1
]
对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种不同的数值算法,可以求解出矩阵指数
eA
e
A
数值运算: E = expm(A)
符号运算:E = expm(sym(A))
矩阵的指数函数: F = expm(A*t)
例4-57 矩阵指数
A=⎡⎣⎢⎢⎢⎢⎢⎢−2001−2001−2−501−5⎤⎦⎥⎥⎥⎥⎥⎥
A
=
[
−
2
1
0
0
−
2
1
0
0
−
2
−
5
1
0
−
5
]
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没有提供对矩阵进行三角函数的现成函数
求解任意非线性矩阵函数的函数调用格式
A1
A
1
= funm(A, ‘fun_name’)
局限性:
采用的算法是矩阵的对角化方法,不能含有重特征根,否则不能得出正确结果
不能求出矩阵的三角函数 sin At
例4-59 矩阵的正弦
A=⎡⎣⎢−301−1−32−1−10⎤⎦⎥
A
=
[
−
3
−
1
−
1
0
−
3
−
1
1
2
0
]
求sin A
A = [-3, -1, -1; 0, -3, -1; 1, 2, 0];
B = funm(A, 'sin')
* 矩阵三角函数的解析解*
根据著名的Euler公式
eja=cosa+jsina
e
j
a
=
cos
a
+
j
sin
a
e−ja=cosa−jsina
e
−
j
a
=
cos
a
−
j
sin
a
推导出
cosa=12(eja+e−ja)
cos
a
=
1
2
(
e
j
a
+
e
−
j
a
)
sina=1j2(eja−e−ja)
sin
a
=
1
j
2
(
e
j
a
−
e
−
j
a
)
此公式可以直接用于矩阵的形式
例4-68 矩阵三角函数的解析解
A=⎡⎣⎢⎢⎢−712−12−4−1−102−60−11−14⎤⎦⎥⎥⎥
A
=
[
−
7
2
0
−
1
1
−
4
2
1
2
−
1
−
6
−
1
−
1
−
1
0
4
]
求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=⎡⎣⎢⎢⎢0000100001000010⎤⎦⎥⎥⎥
H
=
[
0
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
]
幂零矩阵性质
H = diag([1 1 1], 1);
for i = 2: 4
H^i
end
幂零矩阵的高次方是零矩阵
一般矩阵函数的运算
将
mi×mi
m
i
×
m
i
Jordan 块记成
Ji=λiImi+Hmi×mi
J
i
=
λ
i
I
m
i
+
H
m
i
×
m
i
其中,
λi
λ
i
为 Jordan 矩阵的
mi
m
i
重特征值
Hmi
H
m
i
为幂零矩阵,即下式成立
k≥mi,Hkmi≡0
k
≥
m
i
,
H
m
i
k
≡
0
矩阵函数
ψ(Ji)=ψ(λi)Imi+⋯+ψ(mi−1)(λi)(mi−1)!Hmi−1mi
ψ
(
J
i
)
=
ψ
(
λ
i
)
I
m
i
+
⋯
+
ψ
(
m
i
−
1
)
(
λ
i
)
(
m
i
−
1
)
!
H
m
i
m
i
−
1
将任意已知矩阵A进行Jordan分解
A=V⎡⎣⎢⎢⎢⎢⎢J1J2⋱Jm⎤⎦⎥⎥⎥⎥⎥V−1
A
=
V
[
J
1
J
2
⋱
J
m
]
V
−
1
则对该矩阵的任意函数运算符如下
ψ(A)=V⎡⎣⎢⎢⎢⎢⎢ψ(J1)ψ(J2)⋱ψ(Jm)⎤⎦⎥⎥⎥⎥⎥V−1
ψ
(
A
)
=
V
[
ψ
(
J
1
)
ψ
(
J
2
)
⋱
ψ
(
J
m
)
]
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=⎡⎣⎢⎢⎢−712−12−4−1−102−60−11−1−4⎤⎦⎥⎥⎥
A
=
[
−
7
2
0
−
1
1
−
4
2
1
2
−
1
−
6
−
1
−
1
−
1
0
−
4
]
求下列矩阵函数
ψ(A)=eAcosAt
ψ
(
A
)
=
e
A
cos
A
t
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)
矩阵乘方计算 –
Ak
A
k
如何结算
Ak
A
k
?
仿照矩阵任意函数求解,利用Jordan分解与幂零矩阵
如果
A=VJV−1,则Ak=VJkV−1
A
=
V
J
V
−
1
,
则
A
k
=
V
J
k
V
−
1
如果
J=λI+Hm
J
=
λ
I
+
H
m
Jk=λkI+kλk−1Hm+k(k−1)2!λk−2H2m+⋯
J
k
=
λ
k
I
+
k
λ
k
−
1
H
m
+
k
(
k
−
1
)
2
!
λ
k
−
2
H
m
2
+
⋯
例4-71 矩阵乘方计算
A=⎡⎣⎢−301−1−32−1−10⎤⎦⎥
A
=
[
−
3
−
1
−
1
0
−
3
−
1
1
2
0
]
% 乘方计算
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))
矩阵乘方 Ak 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=⎡⎣⎢⎢⎢−712−12−4−1−102−60−11−1−4⎤⎦⎥⎥⎥
A
=
[
−
7
2
0
−
1
1
−
4
2
1
2
−
1
−
6
−
1
−
1
−
1
0
−
4
]
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的平方根与
A12345
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))