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

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

 

上一篇:对于过去所犯的错误,最好的道歉是在将来做正确的事


下一篇:《C语言程序设计:问题与求解方法》——3.2节用计算机求解问题的步骤