线性代数中的数值计算问题

1.AX=b

A\b或者inv(A)*b

一、 特殊矩阵的实现

常见的特殊矩阵有零矩阵、幺矩阵、单位矩阵、三角形矩阵等，这类特殊矩阵在线性代数中具有通用性；还有一类特殊矩阵在专门学科中有用，如有名的希尔伯特(Hilbert)矩阵、范德蒙(Vandermonde) 矩阵等。

1.零矩阵:所有元素值为零的矩阵称为零矩阵。零矩阵可以用zeros函数实现。zeros是MATLAB内部函数，使用格式如下：
zeros(m)：产生m× m阶零矩阵；
zeros(m,n)：产生m ×n阶零矩阵，
                     当m=n时等同于zeros(m)；
zeros(size(A))：产生与矩阵A同样大小的零矩阵。

2.幺矩阵:所有元素值为1的矩阵称为幺矩阵。幺矩阵可以用ones函数实现。它的调用格式与zeros函数一样。
【例1】 试用ones分别建立3×2阶幺矩阵、和与前例矩阵A同样大小的幺矩阵。
用ones(3,2) 建立一个3×2阶幺阵：
ones(3,2) % 一个3×2阶幺阵
ans =1     1
     1     1
     1     1

3.单位矩阵:主对角线的元素值为1、其余元素值为0的矩阵称为单位矩阵。它可以用MATLAB内部函数eye建立，使用格式与zeros相同。
4.数量矩阵:主对角线的元素值为一常数d、其余元素值为0的矩阵称为数量矩阵。显然，当d=1时，即为单位矩阵，故数量矩阵可以用eye(m)*d或eye(m,n)*d建立。

5.对角阵:对角线的元素值为常数、其余元素值为0的矩阵称为对角阵。我们可以通过MATLAB内部函数diag，利用一个向量构成对角阵；或从矩阵中提取某对角线构成一个向量。使用格式为：

diag(V), diag(V,k)
设V为具有m个元素的向量，diag(V)将产生一个m×m阶对角阵，其主对角线的元素值即为向量的元素值；diag(V,k)将产生一个n×n(n=m+|k|，k为一整数)阶对角阵，其第k条对角线的元素值即为向量的元素值。注意：当k＞0，则该对角线位于主对角线的上方第k条；当k＜0，该对角线位于主对角线的下方第|k|条；当k=0，则等同于diag(V)。用diag建立的对角阵必是方阵。

6.上三角阵：使用格式为triu(A)、triu(A,k)
设A为m×n阶矩阵，triu(A)将从矩阵A中提取主对角线之上的上三角部分构成一个m × n阶上三角阵；triu(A,k)将从矩阵A中提取主对角线第|k|条对角线之上的上三角部分构成一个m × n阶上三角阵。注意：这里的k与diag(A,k)的用法类似，当k＞0，则该对角线位于主对角线的上方第k条；当k＜0，该对角线位于主对角线的下方第|k|条；当k=0，则等同于triu (A)

8.下三角阵：使用格式为tril(A)、tril(A,k)
tril的功能是从矩阵A中提取下三角部分构成下三角阵。用法与triu相同。

9.空矩阵
在MATLAB里，把行数、列数为零的矩阵定义为空矩阵。空矩阵在数学意义上讲是空的，但在MATLAB里确是很有用的。例如
A=[0.1 0.2 0.3;0.4 0.5 0.6];
B=find(A>1.0) %返回向量A中符合条件的元素的位置
B = [ ]
这里[ ]是空矩阵的符号，B=find(A>1.0)表示列出矩阵A中值大于1.0的元素的序号。当不能满足括号中的条件时，返回空矩阵。另外，也可以将空矩阵赋给一个变量，如：
B=[ ]
B = [ ]

 二、矩阵的特征值    与特征向量

对于N×N阶方阵A，所谓A的特征值问题是：求数λ和N维非零向量x（通常为复数），使之满足下式：
                         Ax=λx
则称λ为矩阵A的一个特征值（特征根），而非零向量x为矩阵A的特征值λ所对应的特征向量。
对一般的N×N阶方阵A，其特征值通常为复数，若A为实对称矩阵，则A的特征值为实数。

MATLAB提供的内部函数eig可以用来计算特征值与特征向量.eig函数的使用格式有五种,其中常见的有
             E=eig(A),
             [V,D]=eig(A),

(1) E=eig(A)：由eig(A)返回方阵A的N个特征值，构成向量E；
(2) [V,D]=eig(A)：由eig(A)返回方阵A的N个特征值，构成N×N阶对角阵D，其对角线上的N个元素即为相应的特征值，同时将返回相应的特征向量赋予N×N阶方阵V的对应列;

 三、行列式的值

MATLAB提供的内部函数det用来计算矩阵的行列式的值。设矩阵A为一方阵(必须是方阵)，求矩阵A的行列式值的格式为：det(A)。注意：本函数同样能计算通过构造出的稀疏矩阵的行列式的值。

【例6】利用随机函数产生一个三阶方阵A，然后计算方阵之行列式的值。
A=rand(3)
A =
    0.9501    0.4860    0.4565
    0.2311    0.8913    0.0185
    0.6068    0.7621    0.8214
det(A)
ans =
    0.4289

四、   矩阵求逆及其   线性代数方程组求解

1.矩阵的基本性质

① 矩阵的秩：矩阵线性无关的行数与列数称为矩阵的秩。在MATLAB中，求矩阵秩的函数是rank(A)。       ② 矩阵的迹：等于矩阵的对角线元素之和，也等于矩阵的特征值之和。在MATLAB中，求矩阵的迹的函数是trace(A)。
③向量的范数：用来度量矩阵或向量在某种意义下的长度。范数有多种方法定义，其定义不同，范数值也就不同。
(1) norm(V)或norm(V,2)：计算向量V的2—范数。      (2) norm(V,1)：计算向量V的1—范数。      (3) norm(V,inf)：计算向量V的∞—范数。
④ 矩阵的范数：MATLAB提供了求3种矩阵范数的函数，其函数调用格式与求向量的范数的函数完全相同。
⑤  矩阵的条件数：在MATLAB中，计算矩阵A的3种条件数的函数是： (1) cond(A,1)：   计算A的1—范数下的条件数。 (2) cond(A)或cond(A,2)：   计算A的2—范数数下的条件数。 (3) cond(A,inf)：   计算A的 ∞—范数下的条件数。

矩阵条件数

　　矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond(A)=‖A‖·‖A-1‖，对应矩阵的3种范数，相应地可以定义3种条件数。 函数 cond(A,1)、cond(A)或cond(A) 是判断矩阵病态与否的一种度量，条件数越大矩阵越病态。

 

　　条件数事实上表示了矩阵计算对于误差的敏感性。对于线性方程组Ax=b，如果A的条件数大，b的微小改变就能引起解x较大的改变，数值稳定性差。如果A的条件数小，b有微小的改变，x的改变也很微小，数值稳定性好。它也可以表示b不变，而A有微小改变时，x的变化情况。

2 . 矩阵求逆
若方阵A,B满足等式
A*B = B*A = I  (I为单位矩阵)
则称A为B的逆矩阵，或称B为A的逆矩阵。这时A，B都称为可逆矩阵(或非奇异矩阵、或满秩矩阵)，否则称为不可逆矩阵(或奇异矩阵、或降秩矩阵)。

【例7】试用inv函数求方阵A的逆阵A-1赋值给B，且验证A与A-1是互逆的。
A=[1 -1 1;5 -4 3;2 1 1];
B=inv(A)
B =
   -1.4000    0.4000    0.2000
    0.2000   -0.2000    0.4000
    2.6000   -0.6000    0.2000
A*B
ans =
    1.0000    0.0000    0.0000
    0.0000    1.0000    0.0000
    0.0000    0.0000    1.0000

五 、多项式运算及其求根

 MATLAB语言把多项式表达成一个行向量。鉴于MATLAB无零下标，故把多项式的一般形式表达为：

 在MATLAB里，多项式由一个行向量表示，该向量中的元素是按多项式降幂排列的。
              P=[a1  a2  …… an  an+1]
         注意，必须包括具有零系数的项 。

1.多项式求根

命令格式：x=roots(A)。这里A为多项式的系数A(1),A(2),…,A(N),A(N+1)；解得的根赋值给数组X，即X(1),X(2), …,X(N)。
【例9】试用ROOTS函数求多项式x4+8x3-10的根
这是一个4次多项式，它的五个系数依次为：1,8,0,0,-10。下面先产生多项式系数的向量A，然后求根：
A=[1 8 0 0 -10]
A =
     1     8     0     0   -10
x=roots(A)

x =
  -8.0194         
  -0.5075 + 0.9736i
  -0.5075 - 0.9736i
   1.0344

 2.多项式的建立

若已知多项式的全部根，则可以用POLY函数建立起该多项式；也可以用POLY函数求矩阵的特征多项式。POLY函数是一个MATLAB程序，调用它的命令格式是：
                  A=poly(x)
若x为具有N个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋值给向量A。

例:a=[1 2 3;4 5 6;7 8 0];
     p=poly(a)
     p1=poly2str(p,'x')   %显示数学多项式的形式

3.求多项式的值

POLYVAL函数用来求代数多项式的值，调用的命令格式为：
               Y=polyval(A,x)
本命令将POLYVAL函数返回的多项式的值赋值给Y。若x为一数值，则Y也为一数值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。

【例10】以例9的4次多项式、分别取x=1.2和下面的矩阵的2×3个元素为自变量计算该多项式的值。
A=[1 8 0 0 -10];             % 例9的4次多项式系数
x=1.2;                     % 取自变量为一数值
y1=polyval(A,x)
y1 =
  -97.3043
x=[-1 1.2 -1.4;2 -1.8 1.6]      % 给出一个矩阵x
x =
   -1.0000    1.2000   -1.4000

4.多项式的四则运算

(1)多项式加、减
对于次数相同的若干个多项式，可直接对多项式系数向量进行加、减的运算。如果多项式的次数不同，则应该把低次的多项式系数不足的高次项用零补足，使同式中的各多项式具有相同的次数。

(2)多项式乘法
若A、B是由多项式系数组成的向量，则CONV函数将返回这两个多项式的乘积。调用它的命令格式为：
C=conv(A,B)
命令的结果C为一个向量，由它构成一个多项式。

【例11】 a(x)=x2+2x+3; b(x)=4x2+5x+6;
                    c(x) = (x2+2x+3)(4x2+5x+6)

(3)多项式除法
当A、B是由多项式系数组成的向量时，DECONV函数用来对两个多项式作除法运算。调用的命令格式为：
[Q,r]=deconv(A,B)
本命令的结果：多项式A除以多项式B获商多项式赋予Q（也为多项式系数向量）；获余项多项式赋予r（其系数向量的长度与被除多项式相同，通常高次项的系数为0）。
DECONV是CONV的逆函数，即有A=conv(B,Q)+r。

【例12】试用例9的4次多项式与多项式2x2-x+3相除。
A=[1 8 0 0 -10];
B=[2 -1 3];
[P,r]=deconv(A,B)
P =    0.5000    4.2500    1.3750
r =    0         0         0  -11.3750  -14.1250
商多项式P为 0.5x2+4.25x+1.375，
余项多项式r为 -11.375x-14.125。

 (4)多项式的求导
matlab提供了polyder函数多项式的微分。
命令格式：
polyder(p): 求p的微分
polyder(a,b): 求多项式a,b乘积的微分
例：a=[1 2 3 4 5]; poly2str(a,'x')
ans = x^4 + 2 x^3 + 3 x^2 + 4 x + 5
b=polyder(a)
b = 4     6     6     4
poly2str(b,'x')
ans =4 x^3 + 6 x^2 + 6 x + 4