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

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

————————————————————————————————作者：————————————————————————————————日期:

北京航空航天大学

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

第一大题

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

专业: 控制科学与工程

学生姓名：

学号:

教师:

电话:

完成日期: 2015年11月6日

北京航空航天大学

Ｂeｉjing Univeｒsiｔy of Ａeｒｏnaｕtｉcs and Ａstroｎautｉcs

实习题目：

第一题 设有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.求Ａ的与数40

1

5011λλλμ-+=k

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

３.求Ａ的(谱范数)条件数2)A (cond 和行列式d ｅtA 。

说明：

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

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

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

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

一、算法设计方案

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

由于||min ||,501

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

λ其中之一，故可用幂法求出绝对值最大的特征值λ,如果λ=0,则1λ＝λ，否则501λ＝λ。将矩阵A 进行一下平移：

I -A A'λ=

（1）

对'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、求Ａ的（谱范数)条件数2)(A cond 和行列式det Ａ。

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

|

|

)(min max

2λλ=A cond

（3）

其中max λ为按模最大特征值,min λ为按模最小特征值，通过第一问我们求得的

λ和s λ可以很容易求得A 的条件数。

在进行反幂法求解时,要对A 进行LU 分解得到。因L 为单位下三角阵，行列式为1，Ｕ为上三角阵，行列式为主对角线乘积,所以A 的行列式等于U 的行列式,为U 的主对角线的乘积。

二、 算法实现

1、矩阵存储

原矩阵Ａ为一个上、下半带宽都为2的5０1×５01的带状矩阵，由于矩阵中的０元素太多，如果分配一个5０１×501的空间保存矩阵的话会浪费很多空间。因此，为了节省存储量,Ａ的带外元素不给存储，值存储带内元素,如下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 进行D ｏｏl ｉt ｔl ｅ分解，由于A 是带状结构，所以分解出的L 、Ｕ也是带状结构，利用Ｃ矩阵进行Dool ｉt ｔｌｅ分解并求向量k u 的算法如下： （1）作分解A=L Ｕ ﻩ 对于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 的下标都减1。