1.2.14 特殊运算
1．矩阵对角线元素的抽取
函数 diag
格式 X = diag(v,k) %以向量v的元素作为矩阵X的第k条对角线元素，当k=0时，v为X的主对角线；当k>0时，v为上方第k条对角线；当k<0时，v为下方第k条对角线。
X = diag(v) %以v为主对角线元素，其余元素为0构成X。
v = diag(X,k) %抽取X的第k条对角线元素构成向量v。k=0：抽取主对角线元素；k>0：抽取上方第k条对角线元素；k<0抽取下方第k条对角线元素。
v = diag(X) %抽取主对角线元素构成向量v。
例1-46
>> v=[1 2 3];
>> x=diag(v,-1)
x =
0 0 0 0
1 0 0 0
0 2 0 0
0 0 3 0
>> A=[1 2 3;4 5 6;7 8 9]
A =
1 2 3
4 5 6
7 8 9
>> v=diag(A,1)
v =
2
6
2．上三角阵和下三角阵的抽取
函数 tril %取下三角部分
格式 L = tril(X) %抽取X的主对角线的下三角部分构成矩阵L
L = tril(X,k) %抽取X的第k条对角线的下三角部分；k=0为主对角线；k>0为主对角线以上；k<0为主对角线以下。
函数 triu %取上三角部分
格式 U = triu(X) %抽取X的主对角线的上三角部分构成矩阵U
U = triu(X,k) %抽取X的第k条对角线的上三角部分；k=0为主对角线；k>0为主对角线以上；k<0为主对角线以下。
例1-47
>> A=ones(4) %产生4阶全1阵
A =
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
>> L=tril(A,1) %取下三角部分
L =
1 1 0 0
1 1 1 0
1 1 1 1
1 1 1 1
>> U=triu(A,-1) %取上三角部分
U =
1 1 1 1
1 1 1 1
0 1 1 1
0 0 1 1
3．矩阵的变维
矩阵的变维有两种方法，即用“：”和函数“reshape”，前者主要针对2个已知维数矩阵之间的变维操作；而后者是对于一个矩阵的操作。
（1）“：”变维
例1-48
> A=[1 2 3 4 5 6;6 7 8 9 0 1]
A =
1 2 3 4 5 6
6 7 8 9 0 1
>> B=ones(3,4)
B =
1 1 1 1
1 1 1 1
1 1 1 1
>> B(:)=A(:)
B =
1 7 4 0
6 3 9 6
2 8 5 1
（2）Reshape函数变维
格式 B = reshape(A,m,n) %返回以矩阵A的元素构成的m×n矩阵B
B = reshape(A,m,n,p,…) %将矩阵A变维为m×n×p×…
B = reshape(A,[m n p…]) %同上
B = reshape(A,siz) %由siz决定变维的大小，元素个数与A中元素个数
相同。
例1-49 矩阵变维
>> a=[1:12];
>> b=reshape(a,2,6)
b =
1 3 5 7 9 11
2 4 6 8 10 12
4．矩阵的变向
（1）矩阵旋转
函数
格式 B = rot90 (A) %将矩阵A逆时针方向旋转90°
B = rot90 (A,k) %将矩阵A逆时针方向旋转(k×90°)，k可取正负整数。
例1-50
>> A=[1 2 3;4 5 6;7 8 9]
A =
1 2 3
4 5 6
7 8 9
>> Y1=rot90(A),Y2=rot90(A,-1)
Y1 = %逆时针方向旋转
3 6 9
2 5 8
1 4 7
Y2 = %顺时针方向旋转
7 4 1
8 5 2
9 6 3
（2）矩阵的左右翻转
函数 fliplr
格式 B = fliplr(A) %将矩阵A左右翻转
（3）矩阵的上下翻转
函数 flipud
格式 B = flipud(A) %将矩阵A上下翻转
例1-51
>> A=[1 2 3;4 5 6]
A =
1 2 3
4 5 6
>> B1=fliplr(A),B2=flipud(A)
B1 =
3 2 1
6 5 4
B2 =
4 5 6
1 2 3
（4）按指定维数翻转矩阵
函数 flipdim
格式 B = flipdim(A,dim) % flipdim(A,1) = flipud(A)，并且flipdim(A,2)=fliplr(A)。
例1-52
>> A=[1 2 3;4 5 6]
A =
1 2 3
4 5 6
>> B1=flipdim(A,1),B2=flipdim(A,2)
B1 =
4 5 6
1 2 3
B2 =
3 2 1
6 5 4
（5）复制和平铺矩阵
函数 repmat
格式 B = repmat(A,m,n) %将矩阵A复制m×n块，即B由m×n块A平铺而成。
B = repmat(A,[m n]) %与上面一致
B = repmat(A,[m n p…]) %B由m×n×p×…个A块平铺而成
repmat(A,m,n) %当A是一个数a时，该命令产生一个全由a组成的m×n矩阵。
例1-53
>> A=[1 2;5 6]
A =
1 2
5 6
>> B=repmat(A,3,4)
B =
1 2 1 2 1 2 1 2
5 6 5 6 5 6 5 6
1 2 1 2 1 2 1 2
5 6 5 6 5 6 5 6
1 2 1 2 1 2 1 2
5 6 5 6 5 6 5 6
5．矩阵的比较关系
矩阵的比较关系是针对于两个矩阵对应元素的，所以在使用关系运算时，首先应该保证两个矩阵的维数一致或其中一个矩阵为标量。关系运算是对两个矩阵的对应运算进行比较，若关系满足，则将结果矩阵中该位置元素置为1，否则置0。
MATLAB的各种比较关系运算有见表1-1。
表1-1
|
运算符 |
含义 |
运算符 |
含义 |
|
> |
大于关系 |
< |
大于关系 |
|
= = |
等于关系 |
>= |
大于或等于关系 |
|
<= |
小于或等于关系 |
~ = |
不等于关系 |
例1-54
>> A=[1 2 3 4;5 6 7 8];B=[0 2 1 4;0 7 7 2];
>> C1=A==B, C2=A>=B, C3=A~=B
C1 =
0 1 0 1
0 0 1 0
C2 =
1 1 1 1
1 0 1 1
C3 =
1 0 1 0
1 1 0 1
6．矩阵元素的数据变换
对于小数构成的矩阵A来说，如果我们想对它取整数，有以下几种方法：
（1）按-∞方向取整
函数 floor
格式 floor(A) %将A中元素按-∞方向取整，即取不足整数。
（2）按+∞方向取整
函数 ceil
格式 ceil(A) %将A中元素按+∞方向取整，即取过剩整数。
（3）四舍五入取整
函数 round
格式 round (A) %将A中元素按最近的整数取整，即四舍五入取整。
（4）按离0近的方向取整
函数 fix
格式 fix (A) %将A中元素按离0近的方向取整
例1-55
>> A=-1.5+4*rand(3)
A =
2.3005 0.4439 0.3259
-0.5754 2.0652 -1.4260
0.9274 1.5484 1.7856
>> B1=floor(A),B2=ceil(A),B3=round(A),B4=fix(A)
B1 =
2 0 0
-1 2 -2
0 1 1
B2 =
3 1 1
0 3 -1
1 2 2
B3 =
2 0 0
-1 2 -1
1 2 2
B4 =
2 0 0
0 2 -1
0 1 1
（5）矩阵的有理数形式
函数 rat
格式 [n,d]=rat (A) %将A表示为两个整数矩阵相除，即A=n./d。
例1-56 对于上例中的A
>> [n,d]=rat(A)
n =
444 95 131
-225 2059 -472
166 48 1491
d =
193 214 402
391 997 331
179 31 835
（6）矩阵元素的余数
函数 rem
格式 C = rem (A, x) %表示A矩阵除以模数x后的余数。若x=0，则定义rem(A, 0)=NaN，若x≠0，则整数部分由fix(A./x)表示，余数C=A-x.*fix (A./x)。允许模x为小数。
7．矩阵逻辑运算
设矩阵A和B都是m×n矩阵或其中之一为标量，在MATLAB中定义了如下的逻辑运算：
（1）矩阵的与运算
格式 A&B或and(A, B)
说明 A与B对应元素进行与运算，若两个数均非0，则结果元素的值为1，否则为0。
（2）或运算
格式 A|B或or(A, B)
说明 A与B对应元素进行或运算，若两个数均为0，则结果元素的值为0，否则为1。
（3）非运算
格式 ~A或not (A)
说明 若A的元素为0，则结果元素为1，否则为0。
（4）异或运算
格式 xor (A,B)
说明 A与B对应元素进行异或运算，若相应的两个数中一个为0，一个非0，则结果为0，否则为1。
例1-57
>> A=[0 2 3 4;1 3 5 0],B=[1 0 5 3;1 5 0 5]
A =
0 2 3 4
1 3 5 0
B =
1 0 5 3
1 5 0 5
>> C1=A&B,C2=A|B,C3=~A,C4=xor(A,B)
C1 =
0 0 1 1
1 1 0 0
C2 =
1 1 1 1
1 1 1 1
C3 =
1 0 0 0
0 0 0 1
C4 =
1 1 0 0
0 0 1 1
1.2.15 符号矩阵运算
1．符号矩阵的四则运算
Matlab 6.x 抛弃了在4.2版中为符号矩阵设计的复杂函数形式，把符号矩阵的四则运算简化为与数值矩阵完全相同的运算方式，其运算符为：加（＋）、减（－）、乘（×）、除（/、\）等或：符号矩阵的和（symadd）、差（symsub）、乘 (symmul)。
例1-58 >>
>> ；
>>C=B-A
>>D=a\b
则显示：
C=
x-1/x 1-1/(x+1)
x+2-1/(x+2) -1/(x+3)
D=
-6*x-2*x^3-7*x^2 1/2*x^3+x+3/2*x^2
6+2*x^3+10*x^2+14*x -2*x^2-3/2*x-1/2*x^3
2．其他基本运算
符号矩阵的其他一些基本运算包括转置（\'）、行列式（det）、逆（inv）、秩（rank）、幂（^）和指数（exp和expm）等都与数值矩阵相同
3．将数值矩阵转化为符号矩阵
函数 sym
格式 B=sym(A) %将A转化为符号矩阵B
例1-59
>> A=[2/3,sqrt(2),0.222;1.4,1/0.23,log(3)]
A =
0.6667 1.4142 0.2220
1.4000 4.3478 1.0986
>> B=sym(A)
B =
[ 2/3, sqrt(2), 111/500]
[ 7/5, 100/23, 4947709893870346*2^(-52)]
4．符号矩阵的索引与修改
符号矩阵的索引与修改同数值矩阵的索引与修改完全相同，即用矩阵的坐标括号表达式实现。
例1-60 对上例中的矩阵B
>> B(2,3) %矩阵的索引
ans =
4947709893870346*2^(-52)
>> B(2,3)=\'log(7)\' %矩阵的修改
B =
[ 2/3, sqrt(2), 111/500]
[ 7/5, 100/23, log(7)]
5．符号矩阵的简化
符号工具箱中提供了符号矩阵因式分解、展开、合并、简化及通分等符号操作函数。
（1）因式分解
函数 factor
格式 factor(s) %符号表达式s的因式分解函数
说明 S为符号矩阵或符号表达式，常用于多项式的因式分解。
例1-61 将x 9-1分解因式
在Matlab命令窗口键入：
syms x
factor(x^9-1)
则显示：ans =
(x-1)*(x^2+x+1)*(x6+x^3+1)
例1-62 问“入”取何值时，齐次方程组 有非0解？
解：在Matlab编辑器中建立M文件：
syms k
A=[1-k -2 4;2 3-k 1;1 1 1-k];
D=det(A)
factor(D)
其结果显示如下：
D =
-6*k+5*k^2-k^3
ans =
-k*(k-2)*(-3+k)
从而得到：当k=0、k=2或k=3时，原方程组有非0解。
（2）符号矩阵的展开
函数 expand
格式：expand(s) %符号表达式s的展开函数
说明：s为符号矩阵或表达式。常用在多项式的因式分解中，也常用于三角函数，指数函数和对数函数的展开中。
例1-63 将(x+1)3、sin(x+y)展开
在Matlab编辑器中建立M文件：
syms x y
p=expand((x+1)^3)
q=expand(sin(x+y))
则结果显示为
p =
x^3+3*x^2+3*x+1
q =
sin(x)*cos(y)+cos(x)*sin(y)
（3）同类式合并
函数 Collect
格式 Collect(s,v) %将s中的变量v的同幂项系数合并
Collect(s) % s是矩阵或表达式，此命令对由命令findsym函数返回的默认变量进行同类项合并。
（4）符号简化
函数 simple或simplify %寻找符号矩阵或符号表达式的最简型
格式 simple (s) % s是矩阵或表达式
[R,how]=simple (s) %R为返回的最简形，how为简化过程中使用的主要方法。
说明 Simple(s)将表达式s的长度化到最短。若还想让表达式更加精美，可使用函数Pretty。
格式 Pretty(s) %使表达式s更加精美
例1-64 计算行列式 的值。
在Matlab编辑器中建立M文件：
syms a b c d
A=[1 1 1 1;a b c d;a^2 b^2 c^2 d^2;a^4 b^4 c^4 d^4];
d1=det(A)
d2=simple(d1) %化简表达式d1
pretty(d2) %让表达式d2符合人们的书写习惯
则显示结果如下：
d1 =
b*c^2*d^4-b*d^2*c^4-b^2*c*d^4+b^2*d*c^4+b^4*c*d^2-b^4*d*c^2-a*c^2*d^4+a*d^2*c^4+a*b^2*d^4-a*b^2*c^4-a*b^4*d^2+a*b^4*c^2+a^2*c*d^4-a^2*d*c^4-a^2*b*d^4+a^2*b*c^4+a^2*b^4*d-a^2*b^4*c-a^4*c*d^2+a^4*d*c^2+a^4*b*d^2-a^4*b*c^2-a^4*b^2*d+a^4*b^2*c
d2 =
(-d+c)*(b-d)*(b-c)*(-d+a)*(a-c)*(a-b)*(a+c+d+b)
(-d+c)(b-d)(b-c)(-d+a)(a-c)(a-b)(a+c+d+b)
1.2.16 矩阵元素个数的确定
函数 numel
格式 n = numel(a) %计算矩阵A中元素的个数
例1-65
>> A=[1 2 3 4;5 6 7 8];
>> n=numel(A)
n =
8
1.3 矩阵分解
1.3.1 Cholesky分解
函数 chol
格式 R = chol(X) %如果X为n阶对称正定矩阵，则存在一个实的非奇异上三角阵R，满足R\'*R = X；若X非正定，则产生错误信息。
[R,p] = chol(X) %不产生任何错误信息，若X为正定阵，则p=0，R与上相同；若X非正定，则p为正整数，R是有序的上三角阵。
例1-66
>> X=pascal(4) %产生4阶pascal矩阵
X =
1 1 1 1
1 2 3 4
1 3 6 10
1 4 10 20
>> [R,p]=chol(X)
R =
1 1 1 1
0 1 2 3
0 0 1 3
0 0 0 1
p =
0
1.3.2 LU分解
矩阵的三角分解又称LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。
函数 lu
格式 [L,U] = lu(X) %U为上三角阵，L为下三角阵或其变换形式，满足LU=X。
[L,U,P] = lu(X) %U为上三角阵，L为下三角阵，P为单位矩阵的行变换矩阵，满足LU=PX。
例1-67
>> A=[1 2 3;4 5 6;7 8 9];
>> [L,U]=lu(A)
L =
0.1429 1.0000 0
0.5714 0.5000 1.0000
1.0000 0 0
U =
7.0000 8.0000 9.0000
0 0.8571 1.7143
0 0 0.0000
>> [L,U,P]=lu(A)
L =
1.0000 0 0
0.1429 1.0000 0
0.5714 0.5000 1.0000
U =
7.0000 8.0000 9.0000
0 0.8571 1.7143
0 0 0.0000
P =
0 0 1
1 0 0
0 1 0
1.3.3 QR分解
将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。
函数 qr
格式 [Q,R] = qr(A) %求得正交矩阵Q和上三角阵R，Q和R满足A=QR。
[Q,R,E] = qr(A) %求得正交矩阵Q和上三角阵R，E为单位矩阵的变换形式，R的对角线元素按大小降序排列，满足AE=QR。
[Q,R] = qr(A,0) %产生矩阵A的“经济大小”分解
[Q,R,E] = qr(A,0) %E的作用是使得R的对角线元素降序，且Q*R=A(:, E)。
R = qr(A) %稀疏矩阵A的分解，只产生一个上三角阵R，满足R\'*R = A\'*A，这种方法计算A\'*A时减少了内在数字信息的损耗。
[C,R] = qr(A,b) %用于稀疏最小二乘问题：minimize||Ax-b||的两步解：[C,R] = qr(A,b)，x = R\c。
R = qr(A,0) %针对稀疏矩阵A的经济型分解
[C,R] = qr(A,b,0) %针对稀疏最小二乘问题的经济型分解
例1-68
>>A =[ 1 2 3;4 5 6; 7 8 9; 10 11 12];
>>[Q,R] = qr(A)
Q =
-0.0776 -0.8331 0.5444 0.0605
-0.3105 -0.4512 -0.7709 0.3251
-0.5433 -0.0694 -0.0913 -0.8317
-0.7762 0.3124 0.3178 0.4461
R =
-12.8841 -14.5916 -16.2992
0 -1.0413 -2.0826
0 0 0.0000
0 0 0
函数 qrdelete
格式 [Q,R] = qrdelete(Q,R,j) %返回将矩阵A的第j列移去后的新矩阵的qr分解
例1-69
>> A=[-149 -50 -154;537 180 546;-27 -9 -25];
>> [Q,R]=qr(A)
Q =
-0.2671 -0.7088 0.6529
0.9625 -0.1621 0.2176
-0.0484 0.6865 0.7255
R =
557.9418 187.0321 567.8424
0 0.0741 3.4577
0 0 0.1451
>> [Q,R]=qrdelete(Q,R,3) %将A的第3列去掉后进行qr分解。
Q =
-0.2671 -0.7088 0.6529
0.9625 -0.1621 0.2176
-0.0484 0.6865 0.7255
R =
557.9418 187.0321
0 0.0741
0 0
函数 qrinsert
格式 [Q,R] = qrinsert(Q,R,j,x) %在矩阵A中第j列插入向量x后的新矩阵进行qr分解。若j大于A的列数，表示在A的最后插入列x。
例1-70
>> A=[-149 -50 -154;537 180 546;-27 -9 -25];
>> x=[35 10 7]\';
>> [Q,R]=qrinsert(Q,R,4,x)
Q =
-0.2671 -0.7088 0.6529
0.9625 -0.1621 0.2176
-0.0484 0.6865 0.7255
R =
557.9418 187.0321 567.8424 -0.0609
0 0.0741 3.4577 -21.6229
0 0 0.1451 30.1073
1.3.4 Schur分解
函数 schur
格式 T = schur(A) %产生schur矩阵T，即T的主对角线元素为特征值的三角阵。
T = schur(A,flag) %若A有复特征根，则flag=\'complex\'，否则flag=\'real\'。
[U,T] = schur(A,…) %返回正交矩阵U和schur矩阵T，满足A = U*T*U\'。
例1-71
>> H = [ -149 -50 -154; 537 180 546; -27 -9 -25 ];
>> [U,T]=schur(H)
U =
0.3162 -0.6529 0.6882
-0.9487 -0.2176 0.2294
0.0000 0.7255 0.6882
T =
1.0000 -7.1119 -815.8706
0 2.0000 -55.0236
0 0 3.0000
1.3.5 实Schur分解转化成复Schur分解
函数 rsf2csf
格式 [U,T] = rsf2csf (U,T) %将实舒尔形式转化成复舒尔形式
例1-72
>> A=[1 1 1 3;1 2 1 1;1 1 3 1;-2 1 1 4];
>> [u,t]=schur (A)
u =
-0.4916 -0.4900 -0.6331 -0.3428
-0.4980 0.2403 -0.2325 0.8001
-0.6751 0.4288 0.4230 -0.4260
-0.2337 -0.7200 0.6052 0.2466
t =
4.8121 1.1972 -2.2273 -1.0067
0 1.9202 -3.0485 -1.8381
0 0.7129 1.9202 0.2566
0 0 0 1.3474
>> [U,T]=rsf2csf (u,t)
U =
-0.4916 -0.2756 - 0.4411i 0.2133 + 0.5699i -0.3428
-0.4980 -0.1012 + 0.2163i -0.1046 + 0.2093i 0.8001
-0.6751 0.1842 + 0.3860i -0.1867 - 0.3808i -0.4260
-0.2337 0.2635 - 0.6481i 0.3134 - 0.5448i 0.2466
T =
4.8121 -0.9697 + 1.0778i -0.5212 + 2.0051i -1.0067
0 1.9202 + 1.4742i 2.3355 0.1117 + 1.6547i
0 0 1.9202 - 1.4742i 0.8002 + 0.2310i
0 0 0 1.3474
1.3.6 特征值分解
函数 eig
格式 d = eig(A) %求矩阵A的特征值d，以向量形式存放d。
d = eig(A,B) %A、B为方阵，求广义特征值d，以向量形式存放d。
[V,D] = eig(A) %计算A的特征值对角阵D和特征向量V，使AV=VD成立。
[V,D] = eig(A,\'nobalance\') %当矩阵A中有与截断误差数量级相差不远的值时，该指令可能更精确。\'nobalance\'起误差调节作用。
[V,D] = eig(A,B) %计算广义特征值向量阵V和广义特征值阵D，满足AV=BVD。
[V,D] = eig(A,B,flag) % 由flag指定算法计算特征值D和特征向量V，flag的可能值为：\'chol\' 表示对B使用Cholesky分解算法，这里A为对称Hermitian矩阵，B为正定阵。\'qz\' 表示使用QZ算法，这里A、B为非对称或非Hermitian矩阵。
说明 一般特征值问题是求解方程： 解的问题。广义特征值问题是求方程： 解的问题。
1.3.7 奇异值分解
函数 svd
格式 s = svd (X) %返回矩阵X的奇异值向量
[U,S,V] = svd (X) %返回一个与X同大小的对角矩阵S，两个酉矩阵U和V，且满足= U*S*V\'。若A为m×n阵，则U为m×m阵，V为n×n阵。奇异值在S的对角线上，非负且按降序排列。
[U,S,V] = svd (X,0) %得到一个“有效大小”的分解，只计算出矩阵U的前n列，矩阵S的大小为n×n。
例1-73
>> A=[1 2;3 4;5 6;7 8];
>> [U,S,V]=svd(A)
U =
-0.1525 -0.8226 -0.3945 -0.3800
-0.3499 -0.4214 0.2428 0.8007
-0.5474 -0.0201 0.6979 -0.4614
-0.7448 0.3812 -0.5462 0.0407
S =
14.2691 0
0 0.6268
0 0
0 0
V =
-0.6414 0.7672
-0.7672 -0.6414
>> [U,S,V]=svd(A,0)
U =
-0.1525 -0.8226
-0.3499 -0.4214
-0.5474 -0.0201
-0.7448 0.3812
S =
14.2691 0
0 0.6268
V =
-0.6414 0.7672
-0.7672 -0.6414
1.3.8 广义奇异值分解
函数 gsvd
格式 [U,V,X,C,S] = gsvd(A,B) %返回酉矩阵U和V、一个普通方阵X、非负对角矩阵C和S，满足A = U*C*X\'，B = V*S*X\'，C\'*C + S\'*S = I (I为单位矩阵)；A和B的列数必须相同，行数可以不同。
[U,V,X,C,S] = gsvd(A,B,0) %含义与前面相似
sigma = gsvd (A,B) %返回广义奇异值sigma
例1-74
>> A=reshape(1:12,3,4) %产生3行4列矩阵，元素由1，2，…，12构成。
A =
1 4 7 10
2 5 8 11
3 6 9 12
>> B=magic(4) %产生4阶魔方阵
B =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
>> [U,V,X,C,S]=gsvd(A,B)
U =
0.4082 0.7071 0.5774
-0.8165 0.0000 0.5774
0.4082 -0.7071 0.5774
V =
0.2607 -0.7950 -0.5000 0.2236
-0.4029 0.3710 -0.5000 0.6708
-0.5452 -0.0530 -0.5000 -0.6708
0.6874 0.4770 -0.5000 -0.2236
X =
0 -9.4340 -17.0587 3.4641
1.8962 8.7980 -17.0587 8.6603
3.7924 8.1620 -17.0587 13.8564
-5.6885 -7.5260 -17.0587 19.0526
C =
0 0.0000 0 0
0 0 0.0829 0
0 0 0 1.0000
S =
1.0000 0 0 0
0 1.0000 0 0
0 0 0.9966 0
0 0 0 0.0000
