北航数值分析计算实习报告一

航空航天大学

《数值分析》计算实习报告

第一大题

学院：自动化科学与电气工程学院

专业：控制科学与工程

学生姓名：

学号：

教师：

电话：

完成日期： 2015年11月6日

航空航天大学

Beijing University of Aeronautics and Astronautics

实习题目：

第一题 设有501501⨯的实对称矩阵A ，

⎥⎥⎥

⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a

其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1

.0-==⋅⋅⋅=--=c b i e i i a i

i 。矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有

||min ||,501

150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤

1.求1λ，501λ和s λ的值。

2.求A 的与数40

1

5011λλλμ-+=k

k 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。

3.求A 的(谱数)条件数2)A (cond 和行列式detA 。

说明：

1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。

2.选择算法时，应使矩阵A 的所有零元素都不储存。

3.打印以下容： （1）全部源程序；

（2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示12位有效数字。

一、算法设计方案

1、求1λ，501λ和s λ的值。

由于||min ||,501

150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤，可知绝对值最大特征值必为1λ和501

λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则

501λ=λ。将矩阵A 进行一下平移：

I -A A'λ= （1）

对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。

s λ为按模最小特征值，||min ||501

1i i s λλ≤≤=，可对A 使用反幂法求得。

2、求A 的与数40

1

5011λλλμ-+=k

k 最接近的特征值)39,...,2,1(=k k i λ。

计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。因此对A 进行平移变换：

)39,,2,1k -A A k k ==（I μ （2）

对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。 3、求A 的(谱数)条件数2)(A cond 和行列式detA 。

由矩阵A 为非奇异对称矩阵可得：

|

|

)(min max

2λλ=A cond （3）

其中max λ为按模最大特征值，min λ为按模最小特征值，通过第一问我们求得的λ和s λ可以很容易求得A 的条件数。

在进行反幂法求解时，要对A 进行LU 分解得到。因L 为单位下三角阵，行

列式为1，U 为上三角阵，行列式为主对角线乘积，所以A 的行列式等于U 的行列式，为U 的主对角线的乘积。

二、 算法实现

1、矩阵存储

原矩阵A 为一个上、下半带宽都为2的501×501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501×501的空间保存矩阵的话会浪费很多空间。因此，为了节省存储量，A 的带外元素不给存储，值存储带元素，如下C 矩阵所示：

⎥⎥⎥⎥⎥

⎥⎦

⎤⎢⎢⎢

⎢⎢

⎢⎣⎡=00

000050150049932

1c

c

c

c b b b b b a a a a a a b b

b b b

c c c c C

（4） C 是一个5×501的矩阵，相比A 大大节省了存储空间，在数组C 中检索矩阵A 的带元素ij a 的方法是：

j 1,s j -i c a A ++=中的元素的带内元素C ij （5）

2、幂法

幂法迭代公式如下：

⎪⎪⎪⎩⎪

⎪⎪⎨⎧====∈-------k T k k

k k k k k u y Ay u u y u u 111

k 11211k n

0/R βηη任取非零向量 （6）

其中λβ=∞

→k k lim ，不断迭代当εβββ≤--||/||1k k k 时即可认为其满足精度要求，令k βλ=。

在程序中计算1u -=k k Ay 时，根据A 矩阵的特点，简化如下：

⎪⎪⎪⎪

⎩⎪⎪⎪⎪

⎨

⎧⋅⋅⋅=+++++-+-=++=+++=+++=++=)

499,,3()

2()1()()()1()2()()501()500()499()501()

501()500()499()498()500()4()3()2()2(C )1()2()3()2()1(C )1(]

][3[]501][3[]500][3[]2][3[]1][3[i i cy i by i y i C i by i cy i u y C by cy u by y C by cy u cy by y by u cy by y u i （7）

3、反幂法

反幂法迭代公式如下：

⎪⎪⎪⎩⎪

⎪⎪⎨⎧====∈-------k

T k k k k k k k u y y Au u y u u 111

k 11211k n

0/R βηη任取非零向量 （8）

当k 足够大时，k

βλ1

s ≈

。

在求解1-=k k y Au 时，可先对A 进行Doolittle 分解，由于A 是带状结构，所以分解出的L 、U 也是带状结构，利用C 矩阵进行Doolittle 分解并求向量k u 的算法如下： （1）作分解A=LU

对于n ,,2,1k =执行：

)]

,min(,,2,1[/)(:)]

,min(,,1,[:,11

)

,,1max(,1,1,1,11

)

,,1max(,1,1,1,1n r k k k i c c

c

c c n s k k k j c

c

c c k

s k s k r i t k s k t t s t i k s k i k s k i k s j r k t j

s j t t s t k j s j k j s j k +++=-

=++=-

=+---=++-++-++-++----=++-++-++-++-∑∑ （9）

由于C 语言中数组下标是从0开始的，所以在程序中矩阵元素c 的下标都减1。