北航数值分析7

合集下载

北航 数值分析第二次大作业(带双步位移的QR方法)

北航 数值分析第二次大作业(带双步位移的QR方法)

一、算法设计方案:按题目要求,本程序运用带双步位移的QR方法求解给定矩阵的特征值,并对每一实特征值,求解其相应的特征向量。

总体思路:1)初始化矩阵首先需要将需要求解的矩阵输入程序。

为了防止矩阵在后面的计算中被破坏保存A[][]。

2)对给定的矩阵进行拟上三角化为了尽量减少计算量,提高程序的运行效率,在对矩阵进行QR分解之前,先进行拟上三角化。

由于矩阵的QR 分解不改变矩阵的结构,所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。

这里用函数void QuasiTriangularization()来实现,函数形参为double型N维方阵double a[][N]。

3)对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。

用子程序void QR_decomposition()来实现,将Q、R设为形参,然后将计算出来的结果传入Q和R,然后求出RQ乘积。

4)对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛,对QR分解引入双步位移,适当选取位移量,可以避免进行复数运算。

为了进一步减少计算量,在每次进行QR分解之前,先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。

若可以,则得到特征值,同时对矩阵进行降阶处理;若不可以,则进行QR分解。

这里用函数intTwoStepDisplacement_QR()来实现。

这是用来存储计算得到的特征值的二维数组。

考虑到特征值可能为复数,因此将所有特征值均当成复数处理。

此函数中,QR分解部分用子函数void QR_decompositionMk()实现。

这里形参有三个,分别用来传递引入双步位移后的Mk阵,A矩阵,以及当前目标矩阵的维数m。

5)计算特征向量得到特征值后,计算实特征值相应的特征向量。

这里判断所得特征值的虚数部分是否为零。

求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。

因此先初始化矩阵Array,计算(A-λI),再求解方程组。

北航数值分析大作业二(纯原创,高分版)

北航数值分析大作业二(纯原创,高分版)
(R_4 ,I_4 )=( 1.590313458807e+000, 0.000000000000e+000)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令 以及最大迭代步数L;
2.若m≤0,则结束计算,已求出A的全部特征值,判断 或 或m≤2是否成立,成立则转3,否则转4;
3.若 ,则得一个特征值 ,m=m-1,降阶;若 ,则计算矩阵:
的特征值得矩阵A的两个特征值,m=m-2,降阶,转2.;
4.若k≤L,成立则令
k=k+1,转2,否则结束计算,为计算出矩阵A的全部特征值;

北航数值分析实验报告

北航数值分析实验报告

北航‎数值‎分析‎实验‎报告‎‎篇一‎:‎北航‎数值‎分析‎报告‎第一‎大题‎《‎数值‎分析‎》计‎算实‎习报‎告‎第一‎大题‎学‎号:‎D‎Y1‎30‎5‎姓名‎:‎指导‎老师‎:‎一、‎题目‎要求‎已‎知5‎01‎*5‎01‎阶的‎带状‎矩阵‎A,‎其特‎征值‎满足‎?1‎?‎2‎..‎.‎?5‎01‎。

试‎求:‎1‎、?‎1,‎?5‎01‎和?‎s的‎值;‎‎2、‎A的‎与数‎?k‎??‎1?‎k‎?5‎01‎??‎1‎40‎最‎接近‎的特‎征值‎?i‎k(‎k=‎1,‎2,‎..‎.,‎39‎);‎‎3、‎A的‎(谱‎范数‎)条‎件数‎c n‎d(‎A)‎2和‎行列‎式d‎e t‎A。

‎‎二、‎算法‎设计‎方案‎题‎目所‎给的‎矩阵‎阶数‎过大‎,必‎须经‎过去‎零压‎缩后‎进行‎存储‎和运‎算,‎本算‎法中‎压缩‎后的‎矩阵‎A1‎如下‎所示‎。

‎?0‎?0‎?A‎1?‎?a‎1‎??‎b?‎?c‎0‎b a‎2b‎c‎c b‎b c‎.‎..‎..‎..‎..‎..‎.‎c b‎b c‎c‎b a‎50‎0b‎0‎a ‎3.‎..‎a4‎99‎c‎?‎b?‎?a‎50‎1?‎?‎0?‎0?‎?‎由矩‎阵A‎的特‎征值‎满足‎的条‎件可‎知‎?1‎与?‎50‎1之‎间必‎有一‎个最‎大,‎则采‎用幂‎法求‎出的‎一‎个特‎征值‎必为‎其中‎的一‎个:‎当‎所求‎得的‎特征‎值为‎正数‎,则‎为?‎50‎1;‎否则‎为?‎1。

‎在求‎得?‎1与‎?‎50‎1其‎中的‎一个‎后,‎采用‎带位‎移的‎幂法‎则可‎求出‎它们‎中的‎另一‎个,‎且位‎移量‎即为‎先求‎出的‎特‎征值‎的值‎。

用‎反幂‎法求‎得的‎特征‎值必‎为?‎s。

‎由条‎件数‎的性‎质可‎得,‎c n‎d(‎A)‎2为‎模最‎大的‎特征‎值与‎模最‎小的‎特征‎值之‎比的‎模,‎因此‎,求‎出?‎1,‎?5‎01‎和?‎s的‎值后‎,则‎可以‎求得‎c n‎d(‎A)‎2。

北航研究生数值分析上机作业 三 (报告+所有程序大全)

北航研究生数值分析上机作业 三 (报告+所有程序大全)

数值分析上机作业3——求解非线性方程组以及二元函数的插值拟合1. 算法设计对于全部的插值节点(,),0,1,...,10,0,1,...,20i j x y i j ==,带入非线性方程组中,用Newton 迭代法解非线性方程组,得到(,),0,1,...,10,0,1,...,20i j t u i j ==。

对(,)i j t u ,在二维数表中进行插值,采用分片双二次插值法。

插值过程中,先选择分片区域的中心节点,在数表中的列记为(0:5)tt ,行记为(0:5)uu ,中心节点记为(,)a b ,生成向量_(0:2)t temp ,_(0)(())((1))/(((1)())((1)(1)))i i t temp t tt a t tt a tt a tt a tt a tt a =--+----+, _(1)((1))((1))/((()(1))(()(1)))i i t temp t tt a t tt a tt a tt a tt a tt a =---+---+, _(2)((1))(())/(((1)(1))((1)()))i i t temp t tt a t tt a tt a tt a tt a tt a =---+--+-,同理,生成向量_(0:2)u temp ,_(0)(())((1))/(((1)())((1)(1)))_(1)((1))((1))/((()(1))(()(1)))_(2)((1))(())/(((1)(1))((1)())j j j j j j u temp u uu a u uu a uu a uu a uu a uu a u temp u uu a u uu a uu a uu a uu a uu a u temp u uu a u uu a uu a uu a uu a uu a =--+----+=---+---+=---+--+-)记数表中以分片区域中心节点为中心的3×3的矩阵为T , 对于(,)i j t u 插值结果为(_)()(_)T t temp T u temp 。

北航研究生数值分析试题

北航研究生数值分析试题

∗⎞ ⎟的 A1 ⎠
矩阵。
三、(12 分)试用高斯列主元素法求解线性方程组
⎡ 1 3 −2 −4 ⎤ ⎡ x1 ⎤ ⎡3 ⎤ ⎢ 2 6 −7 −10 ⎥ ⎢ x ⎥ ⎢ −2 ⎥ ⎢ ⎥⎢ 2⎥ = ⎢ ⎥ ⎢ −1 −1 5 9 ⎥ ⎢ x3 ⎥ ⎢14 ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ x4 ⎦ ⎥ ⎣ −6 ⎦ ⎣ −3 −5 0 15 ⎦ ⎣ 四、(12 分)利用矩阵 A 的三角分解 A = LU 求解下列方程组 ⎛ 1 2 1 ⎞ ⎛ x1 ⎞ ⎛ 0 ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ 2 2 3 ⎟ ⎜ x2 ⎟ = ⎜ 3 ⎟ ⎜ −1 −3 0 ⎟ ⎜ x ⎟ ⎜ 2 ⎟ ⎝ ⎠⎝ 3 ⎠ ⎝ ⎠
第一章
1、近似数 x = 0.231 关于真值 x = 0.229 有( (1)1;(2)2;(3)3;(4)4。

绪论
一、选择题(四个选项中仅有一项符合题目要求,每小题 3 分,共计 15 分) )位有效数字。
2、取 3 ≈ 1.732 计算 x = ( 3 − 1) ,下列方法中哪种最好?(
4

Ax
∞和
A ∞ 的值分别为(

3
(1) 8 , 8 ;
(2) 8 , 7 ;
(3) 8 , 6 ;
(4) 7 , 7 。
5 、若解线性代数方程组的 Gauss 部分选主元方法第二步得到的系数矩阵的第三列向量为
(2
6 3 2 −5 4 2 ) ,则第三步主行是(
T
) (4) 第 6 行。
(1) 第 2 行;
1 − cos x , sin x
x ≠ 0且 x << 1 ;
(2)
1 1− x , − 1+ 2x 1+ x

北航2010-2015年研究生数值分析报告期末模拟试卷与真题

北航2010-2015年研究生数值分析报告期末模拟试卷与真题

北航2010-2015年研究生数值分析报告期末模拟试卷与真题数值分析模拟卷A一、填空(共30分,每空3分)1 设-=1511A ,则A 的谱半径=)(a ρ______,A 的条件数)(1A cond =________. 2 设 ,2,1,0,,53)(2==+=k kh x x x f k ,则],,[21++n n n x x x f =________, ],,[321+++n n n n x x x x f ,=________.3 设≤≤-++≤≤+=21,1210,)(2323x cx bx x x x x x S ,是以0,1,2为节点的三次样条函数,则b=________,c=________.4 设∞=0)]([k k x q 是区间[0,1]上权函数为x x =)(ρ的最高项系数为1的正交多项式族,其中1)(0=x q ,则?=10)(dx x xq k ________,=)(2x q ________.5 设=11001a a a a A ,当∈a ________时,必有分解式,其中L 为下三角阵,当其对角线元素)3,2,1(=i L ii 满足条件________时,这种分解是唯一的.二、(14分)设49,1,41,)(21023====x x x x x f , (1)试求)(x f 在]49,41[上的三次Hermite 插值多项式)(x H 使满足2,1,0),()(==i x f x H i i ,)()(11x f x H '='.(2)写出余项)()()(x H x f x R -=的表达式.三、(14分)设有解方程0cos 2312=+-x x 的迭代公式为n n x x cos 3241+=+,(1)证明R x ∈?0均有?∞→=x x n x lim (?x 为方程的根);(2)取40=x ,用此迭代法求方程根的近似值,误差不超过,列出各次迭代值;(3)此迭代的收敛阶是多少?证明你的结论.四、(16分) 试确定常数A ,B ,C 和,使得数值积分公式有尽可能高的代数精度. 试问所得的数值积分公式代数精度是多少?它是否为Gauss 型的?五、(15分)设有常微分方程的初值问题=='00)(),(y x y y x f y ,试用Taylor 展开原理构造形如)()(11011--++++=n n n n n f f h y y y ββα的方法,使其具有二阶精度,并推导其局部截断误差主项.六、(15分)已知方程组b Ax =,其中= ??=21,13.021b A ,(1)试讨论用Jacobi 迭代法和Gauss-Seidel 迭代法求解此方程组的收敛性.(2)若有迭代公式)()()()1(b Ax a x x k k k ++=+,试确定一个的取值围,在这个围任取一个值均能使该迭代公式收敛.七、(8分)方程组,其中,A 是对称的且非奇异.设A 有误差,则原方程组变化为,其中为解的误差向量,试证明 .其中1λ和2λ分别为A 的按模最大和最小的特征值.数值分析模拟卷B填空题(每空2分,共30分)1. 近似数231.0=*x 关于真值229.0=x 有____________位有效数字;2. 设)(x f 可微,求方程)(x f x =根的牛顿迭代格式是_______________________________________________;3. 对1)(3++=x x x f ,差商=]3,2,1,0[f _________________;=]4,3,2,1,0[f ________;4. 已知???? ??-='-=1223,)3,2(A x ,则=∞||||Ax ________________,=)(1A Cond ______________________ ;5. 用二分法求方程01)(3=-+=x x x f 在区间[0,1]的根,进行一步后根所在区间为_________,进行二步后根所在区间为_________________;6. 求解线性方程组=+=+04511532121x x x x 的高斯—赛德尔迭代格式为_______________________________________;该迭代格式迭代矩阵的谱半径=)(G ρ_______________;7. 为使两点数值求积公式:?-+≈111100)()()(x f x f dx x f ωω具有最高的代数精确度,其求积节点应为=0x _____ , =1x _____,==10ωω__________.8. 求积公式)]2()1([23)(30f f dx x f +≈?是否是插值型的__________,其代数精度为___________。

北航数值分析-lec7-幂法和反幂法

北航数值分析-lec7-幂法和反幂法
线性方程组求解
迭代收敛性
反幂法在求解特征值问题中的应用
特征值问题
反幂法主要用于求解矩阵的特征值和特征向量问题。通过迭代过程,反幂法能够找到矩阵的所有特征 值和对应的特征向量。
数值稳定性
反幂法在求解特征值问题时,需要关注数值稳定性问题。由于计算机浮点运算的误差累积,反幂法可 能会产生数值不稳定的解。因此,需要采取适当的策略来提高数值稳定性。
误差分析比较
幂法
由于幂法是通过连续的矩阵乘法来计算矩阵的幂,因此误差会随着计算的次数逐渐 累积。为了减小误差,需要选择合适的计算精度和减少计算次数。
反幂法
反幂法是通过求解线性方程组来计算矩阵的逆和行列式,因此误差主要来自于线性 方程组的求解精度。为了减小误差,需要选择合适的求解方法和提高求解精度。
202X
北航数值分析-lec7-幂法 和反幂法
单击此处添加副标题内容
汇报人姓名 汇报日期
目 录幂法介绍Fra bibliotek反幂法介绍
幂法和反幂法的比较
幂法和反幂法的实现细节
幂法和反幂法的实际应用案例
单击此处输入你的正文,文字是
您思想的提炼,请尽量言简意赅
的阐述观点
contents
单击此处输入你的正文,文字是 您思想的提炼,请尽量言简意赅 的阐述观点
反幂法的实现细节
反幂法是一种迭代算法,用 于求解线性方程组的近似逆。
反幂法的收敛速度取决于矩阵的谱 半径,如果矩阵的谱半径较小,则 反幂法收敛速度较快。
ABCD
反幂法的实现步骤包括:选择初始 矩阵、计算迭代矩阵、更新解矩阵 和判断收敛性。
在实际应用中,反幂法通常用于 求解大规模稀疏线性系统的预处 理和后处理问题。
01

北京航空航天大学数值分析课程知识点总结

北京航空航天大学数值分析课程知识点总结

1 ,其中 1 和 n 分别是矩阵 A 的 n
2.4 迭代法
2.4.1 迭代法的一般形式及其收敛性
x ( k 1) Gx ( k ) d (k 0,1,...)
定义 设 n n 矩阵 G 的特征值是 1 , 2 ,..., n ,称 (G ) max | i | 为矩阵 G 的谱半径。
n T
x 1 xi
i 1
n
x2 x
则 1 , 2 和 都是向量范数。 定理 1.2 设


x
i 1 1 i n
n
2 i
max xi


是 R 上的任意两种向量范数,则存在与向量 x 无关的常数 m 和
n
M(0<m<M),使下列关系式成立
m x
1.3.2 矩阵范数
~
若 f '(a ) 0 且 | f ''( a ) | / | f '( a ) | 不很大,则有误差估计
e(u ) f '(a )e(a )
~
(u ) f '(a) (a)
~

若 f '(a ) f ''(a ) ... f
( k 1)
(a ) 0, f
(k )
... ... ... ... ln ,n 1
为节省空间,用 C(m,n)存储 A 的带内元素,其中 m=r+s+1,并且 aij ci j s 1, j 。 2.2.5 拟三对角线性方程组的求解方法
a1 d 2 A cn p1 d 2 r1
e xa e ,称 er 为近似值 a 的相对误差。由于 x 未知,实际上总把 作为 a 的 x x a e xa , 相对误差一般用百分比表示。er 的上界, 即r a a |a|

北航《数值分析》习题

北航《数值分析》习题

北航《数值分析》习题习题一1. 下列各数都是经过四舍五入得到的近似数,试指出它们有几位有效数字以及它们的绝对误差限、相对误差限。

(1);(2);(3);(4);(5);(6);(7);1. (1)5,,;(2)2,,;(3)4,,;(4)5,,;(5)1,,;(6)2,,;(7)6,,2. 为使下列各数的近似值的相对误差限不超过,问各近似值分别应取几位有效数字?2. ;;3. 设均为第1题所给数据,估计下列各近似数的误差限。

(1);(2);(3)3. (1);(2);(3)4. 计算,取,利用下列等价表达式计算,哪一个的结果最好?为什么?(1);(2);(3)(4)4. 第(3)个结果最好。

5. 序列满足递推关系式若(三位有效数字),计算时误差有多大?这个计算过程稳定吗?5. 不稳定。

从计算到时,误差约为6. 求方程的两个根,使其至少具有四位有效数字(要求利用。

6. ,7. 某生产部门生产的一件产品需用七个零件,而这七个零件的质量取决于零件参数的标定值,它们的参数允许有一定的误差:若每一零件的标定值取做区间中点,在生产过程中每一零件的参数都有可能产生误差。

由此将零件分成不同的等级:A,B,C三等,等级由标定值的相对误差限表示,A等为1%,B等为5%,C等为10%。

试确定三个等级的零件分别满足的区间。

8. 将一个八位二进制数(10111101)2转换成十进制数时,可以用公式:(1)用多项式求值的秦九韶方法求C的值;(2)写出将任意一个八位二进制数(b1b2b3b4b5b6b7b8)2转化为十进制数的算法。

9. 利用等式变换使下列表达式的计算结果比较精确。

(1);(2)(3);(4)9. (1);(2);(3);(4)10. 设,求证:(1)(2)利用(1)中的公式正向递推计算时误差增大;反向递推时误差函数减小。

习题二1. 判断下列方程有几个实根,并求出其隔根区间。

(1);(2)(3);(4)1. (1),,;(2);(3),,;(4)为根。

北航研究生数值分析编程大作业1

北航研究生数值分析编程大作业1

数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2,设置矩阵[][]5011++s r C ,在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是:j s j i ij c a ,1++-=。

这样所需要的存储单元数大大减少,从而极大提高了运算效率。

2、利用幂法求出5011λλ,幂法迭代格式:011111111nT k k k k k k kk T k k k u R u u y u u Ay y u ηηβ--------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时,迭代终止。

首先对于矩阵A 利用幂法迭代求出一个λ,然后求出矩阵B ,其中I A B λ-=(I 为单位矩阵),对矩阵B 进行幂法迭代,求出λ',之后令λλλ+'='',比较的大小与λλ'',大者为501λ,小者为1λ。

3、利用反幂法求出ik s λλ,反幂法迭代格式:011111111nTk k k k k k kk T k k k u R u u y u Au y y u ηηβ--------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时,迭代终止,1s k λβ=。

每迭代一次都要求解一次线性方程组1-=k k y Au ,求解过程为:(1)作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1(2)求解y Ux b Ly ==,(数组b 先是存放原方程组右端向量,后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法,直接可以求得矩阵按模最小的特征值s λ。

北航 数值分析第二次大作业(带双步位移的QR方法)

北航 数值分析第二次大作业(带双步位移的QR方法)

一、算法设计方案:按题目要求,本程序运用带双步位移的QR方法求解给定矩阵的特征值,并对每一实特征值,求解其相应的特征向量。

总体思路:1)初始化矩阵首先需要将需要求解的矩阵输入程序。

为了防止矩阵在后面的计算中被破坏保存A[][]。

2)对给定的矩阵进行拟上三角化为了尽量减少计算量,提高程序的运行效率,在对矩阵进行QR分解之前,先进行拟上三角化。

由于矩阵的QR 分解不改变矩阵的结构,所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。

这里用函数void QuasiTriangularization()来实现,函数形参为double型N维方阵double a[][N]。

3)对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。

用子程序void QR_decomposition()来实现,将Q、R设为形参,然后将计算出来的结果传入Q和R,然后求出RQ乘积。

4)对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛,对QR分解引入双步位移,适当选取位移量,可以避免进行复数运算。

为了进一步减少计算量,在每次进行QR分解之前,先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。

若可以,则得到特征值,同时对矩阵进行降阶处理;若不可以,则进行QR分解。

这里用函数intTwoStepDisplacement_QR()来实现。

这是用来存储计算得到的特征值的二维数组。

考虑到特征值可能为复数,因此将所有特征值均当成复数处理。

此函数中,QR分解部分用子函数void QR_decompositionMk()实现。

这里形参有三个,分别用来传递引入双步位移后的Mk阵,A矩阵,以及当前目标矩阵的维数m。

5)计算特征向量得到特征值后,计算实特征值相应的特征向量。

这里判断所得特征值的虚数部分是否为零。

求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。

因此先初始化矩阵Array,计算(A-λI),再求解方程组。

北航数值分析大作业第二题

北航数值分析大作业第二题

数值分析第二次大作业史立峰SY1505327一、 方案(1)利用循环结构将sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值,得到需要变换的矩阵A ;(2)然后,对矩阵A 利用Householder 矩阵进行相似变换,把A 化为上三角矩阵A(n-1)。

对A 拟上三角化,得到拟上三角矩阵A(n-1),具体算法如下:记A(1)=A ,并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)(ΛΛ+==。

对于2,,2,1-=n r Λ执行 1. 若()n r r i a r ir,,3,2)(Λ++=全为零,则令A(r+1) =A(r),转5;否则转2。

2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,12r rr r r r a c c h +-=3. 令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0ΛΛ。

4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。

(3)使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值,也是A 的全部特征值,具体算法如下:1. 给定精度水平0>ε和迭代最大次数L 。

2. 记n n ij n a A A ⨯-==][)1()1()1(,令n m k ==,1。

3. 如果ε≤-)(1,k m m a ,则得到A 的一个特征值)(k mm a ,置1:-=m m (降阶),转4;否则转5。

北航数值分析第一次大作业

北航数值分析第一次大作业

一、算法的设计方案:(一)各所求值得计算方法1、最大特征值λ501,最小特征值λ1,按模最小特征值λs的计算方法首先使用一次幂法运算可以得到矩阵的按模最大的特征值λ,λ必为矩阵A的最大或最小特征值,先不做判断。

对原矩阵A进行一次移项,即(A-λI),在进行一次幂法运算,可以得到另一个按模最大特征值λ0。

比较λ和λ的大小,较大的即为λ501,较小的即为λ1。

对矩阵A进行一次反幂法运算,即可得到按模最小特征值λs。

2、A与μk 值最接近的特征值λik的计算方法首先计算出k所对应的μk 值,对原矩阵A进行一次移项,即(A-μkI),得到一个新的矩阵,对新矩阵进行一次反幂法运算,即可得到一个按模最小特征值λi 。

则原矩阵A与μk值最接近的特征值λik=λi+μk。

3、A的(谱范数)条件数cond(A)2的计算方法其中错误!未找到引用源。

矩阵A的按模最大和按模最小特征值。

(二)程序编写思路。

由于算法要求A的零元素不存储,矩阵A本身为带状矩阵,所以本题的赋值,LU分解,反幂法运算过程中,均应采用Doolittle分解法求解带状线性方程组的算法思路。

幂法、反幂法和LU分解均是多次使用,应编写子程序进行反复调用。

二、源程序:#include<stdio.h>#include<iostream>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip> /*头文件*//*定义全局变量*/#define N 502 /*取N为502,可实现从1到501的存储,省去角标变换的麻烦*/ #define epsilon 1.0e-12 /*定义精度*/#define r 2 /*r,s为带状矩阵的半带宽,本题所给矩阵二者都是2*/ #define s 2double c[6][N]; /*定义矩阵存储压缩后的带状矩阵*/double fuzhi(); /*赋值函数*/void LUfenjie(); /*LU分解程序*/int max(int a,int b); /*求两个数字中较大值*/int min(int a,int b); /*求两个数字中较小值*/double mifa(); /*幂法计算程序*/double fanmifa(); /*反幂法计算程序*/double fuzhi() /*赋值程序,按行赋值,行从1到5,列从1到501,存储空间的第一行第一列不使用,角标可以与矩阵一一对应,方便书写程序*/{int i,j;i=1;for(j=3;j<N;j++){c[i][j]=-0.064;}i=2;for(j=2;j<N;j++){c[i][j]=0.16;}i=3;for(j=1;j<N;j++){c[i][j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);}i=4;for(j=1;j<N-1;j++){c[i][j]=0.16;}i=5;for(j=1;j<N-2;j++){c[i][j]=-0.064;}return(c[i][j]);}int max(int a,int b){ return((a>b)?a:b);}int min(int a,int b){ return((a<b)?a:b);}void LUfenjie() /*LU分解程序,采用的是带状矩阵压缩存储后的LU分解法*/{double temp;int i,j,k,t;for(k=1;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++){temp=0;for(t=max(1,max(k-r,j-s));t<=(k-1);t++){temp=temp+c[k-t+s+1][t]*c[t-j+s+1][j];}c[k-j+s+1][j]=c[k-j+s+1][j]-temp;}for(i=k+1;i<=min(k+r,N-1);i++){temp=0;for(t=max(1,max(i-r,k-s));t<=(k-1);t++){temp=temp+c[i-t+s+1][t]*c[t-k+s+1][k];}c[i-k+s+1][k]=(c[i-k+s+1][k]-temp)/c[s+1][k];}}}double mifa() /*幂法计算程序*/ {double u0[N],u1[N];double temp,Lu,beta=0,beta0;int i,j;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;}for(i=1;i<N;i++){temp=0;for(j=max(i-1,1);j<=min(i+2,N-1);j++){temp=temp+c[i-j+s+1][j]*u1[j]; }u0[i]=temp;} //新的u0temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}double fanmifa() /*反幂法计算程序*/{double u0[N],u1[N],u2[N];double temp,Lu,beta=0,beta0;int i,t;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;u2[i]=u1[i];}fuzhi();LUfenjie();/*带状矩阵压缩存储并进行LU分解后,求解线性方程组得到迭代向量u k,即程序中的u0*/for(i=2;i<N;i++){ temp=0;for(t=max(1,i-r);t<=(i-1);t++){temp=temp+c[i-t+s+1][t]*u2[t];}u2[i]=u2[i]-temp;}u0[N-1]=u2[N-1]/c[s+1][N-1];for(i=N-2;i>=1;i--){ temp=0;for(t=i+1;t<=min(i+s,N-1);t++){temp=temp+c[i-t+s+1][t]*u0[t];}u0[i]=(u2[i]-temp)/c[s+1][i];}temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;beta=1/beta; /*beta即为所求特征值,可直接返回*/}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}void main(){double u[40]; /*定义数组,存放k值运算得到的μk值*/double lambda1,lambda501,lambdak,a,b,d,cond,det;int i,j,k;fuzhi();a=mifa(); /*幂法计算按模最大值*/fuzhi();d=fanmifa(); /*反幂法计算按模最小值*/fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-a;}b=mifa()+a; /*移项后幂法计算按模最大值*/if(a>b) /*比较两个按模最大值大小,并相应输出最大特征值λ501和最小特征值λ1*/ {lambda1=b;lambda501=a;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}else{lambda1=a;lambda501=b;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}printf("矩阵A按模最小特征值lambdas=%13.11e\n",d); /*输出按模最小特征值λs*/for(k=1;k<40;k++) /*对每一个进行移项反幂法运算,求出最接近μk的特征值并输出*/ {u[k]=(lambda501-lambda1)*k/40+lambda1;fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-u[k];}lambdak=fanmifa()+u[k];i=k;printf("矩阵A最接近uk特征值lambdak%d=%13.11e\n",i,lambdak);}cond=fabs(a/d);printf("A的条件数=%13.11e\n",cond); /*计算A条件数并输出*/fuzhi(); /*计算A的行列式值并输出*/LUfenjie();det=1;for(i=1;i<N;i++){det=det*c[3][i];}printf("行列式的值detA=%13.11e\n",det);}三、程序的运行结果:四、初始向量的选取对计算结果的影响:(一)选取形式不变,数值变换1、取u0为[0.5,0.5………..0.5],运行结果如下:2、取u0为[50,50………..50],运行结果如下:从运行结果来看,此类初始向量的选取对结果不会产生影响,即使选成0,结果也不变化。

北航数值分析作业第一题

北航数值分析作业第一题

数值分析作业第一题一、 算法设计方案利用带状Dollittle 分解,将A[501][501]转存到数组C[5][501],以节省存储空间1、计算λ1和λ501首先使用幂法求出矩阵的按模最大的特征值λ0:如果λ0>0,则其必为按模最大值,因此λ501=λ0,然后采用原点平移法,平移量为λ501,使用幂法迭代求出矩阵A -λ501I 的按模最大的特征值,其特征值按从小到大排列应为λ1-λ501、λ2-λ501、……、0。

因此A-λ501I 的按模最大的特征值应为λ1-λ501。

又因为λ501的值已求得,由此可直接求出λ1。

2、计算λSλS 为矩阵A 按模最小的特征值,可以通过反幂法直接求出。

3、计算λikλik 是对矩阵A 进行λik 平移后,再用反幂法求出按模最小的特征值λmin ,λik =λik +λmin 。

4、计算矩阵A 的条件数计算cond (A )2和行列式det(A)矩阵A 的条件数为n12cond λλ)( A ,其中λ1和λn 分别是矩阵A 的模最大和最小特征值,直接利用上面求得的结果直接计算。

矩阵A 的行列式可先对矩阵A 进行LU 分解后,det(A)等于U 所有对角线上元素的乘积。

二、源程序:#include<math.h>#include<stdio.h>#include<stdlib.h>#include<iostream.h>#define s 2#define r 2int Max(int v1,int v2);int Min(int v1,int v2);int maxt(int v1,int v2,int v3);void storage(double C[5][501],double b,double c);double mifa(double C[5][501]);void LU(double C[5][501]);double fmifa(double C[5][501]);int Max(int v1,int v2) //求两个数的最大值{ return((v1>v2)?v1:v2);}int Min(int v1,int v2) //求两个数最小值{ return ((v1<v2)?v1:v2);}int maxt(int v1,int v2,int v3) //求三个数最大值{ int t;if(v1>v2) t=v1;else t=v2;if(t<v3) t=v3;return(t);}/***将矩阵值转存在一个数组里,以节省存储空间***/void storage(double C[5][501],double b,double c){ int i=0,j=0;C[i][j]=0,C[i][j+1]=0;for(j=2;j<=500;j++)C[i][j]=c;i++;j=0;C[i][j]=0;for(j=1;j<=500;j++)C[i][j]=b;i++;for(j=0;j<=500;j++)C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1));i++;for(j=0;j<=499;j++)C[i][j]=b;C[i][j]=0;i++;for(j=0;j<=498;j++)C[i][j]=c;C[i][j]=0,C[i][j+1]=0;}//用于求解最大的特征值,幂法double mifa(double C[5][501]){ int m=0,i,j;double b2,b1=0,sum;double u[501],y[501];for (i=0;i<501;i++){ u[i] = 1.0;}do{ sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){ u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+C[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);return b2;}/*****行列式LU分解*****/void LU(double C[5][501]){ double sum;int k,i,j;for(k=1;k<=501;k++){ for(j=k;j<=Min(k+s,501);j++){ sum=0;for(i=maxt(1,k-r,j-s);i<=k-1;i++)sum+=C[k-i+s][i-1]*C[i-j+s][j-1];C[k-j+s][j-1]-=sum;}for(j=k+1;j<=Min(k+r,501);j++){ sum=0;for(i=maxt(1,j-r,k-s);i<=k-1;i++)sum+=C[j-i+s][i-1]*C[i-k+s][k-1];C[j-k+s][k-1]=(C[j-k+s][k-1]-sum)/C[s][k-1];}}}/***带状DOOLITE分解,并且求解出方程组的解***/void solve(double C[5][501],double x[501],double b[501]){ int i,j,k,t;double B[5][501],c[501];for(i=0;i<=4;i++){ for(j=0;j<=500;j++)B[i][j]=C[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){ for(j=k;j<=Min(k+s,500);j++){ for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){ for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k];B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){ x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//用于求解模最大的特征值,反幂法double fmifa(double C[5][501]){ int m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){ [i] = 1.0;}do{ if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);solve(C,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=1.0e-12);return 1/b2;}/***主程序***/void main(){ double b=0.16,c=-0.064,det=1.0;int i;double C[5][501],cond;storage(C,b,c); //进行C的赋值cout.precision(12); //定义输出精度double k1=mifa(C); //利用幂法计算矩阵的最大特征值和最小特征值if(k1<0)printf("λ1=%.12e\n",k1);else if(k1>=0)printf("λ501=%.12e\n",k1);for(i=0;i<501;i++)C[2][i]=C[2][i]-k1;double k2=mifa(C)+k1;if(k2<0)printf("λ1=%.12e\n",k2);else if(k2>=0)printf("λ501=%.12e\n",k2);storage(C,b,c);double k3=fmifa(C); //利用反幂法计算矩阵A的按模最小特征值printf("λs=%.12e\n",k3);storage(C,b,c); //计算最接近特征值double u[39]={0};for(i=0;i<39;i++){ u[i]=k1+(i+1)*(k2-k1)/40;C[2][i]=C[2][i]-u[i];u[i]=fmifa(C)+u[i];printf("与数u%d 最接近的特征值λ%d: %.12e\n",i+1,i+1,u[i]);}if(k1>0) //计算矩阵A的条件数,取2范数cond=fabs(k1/k3);else if(k1<0)cond=fabs(k2/k3);storage(C,b,c);LU(C); //利用LU分解计算矩阵A的行列式for(i=0;i<501;i++)det*=C[2][i];printf("\ncond(A)=%.12e\n",cond);printf("\ndet(A)=%.12e\n",det);}三、计算结果:四、结果分析迭代初始向量的选择对果有一定的影响,选择不同的初始向量可能会得到不同阶的特征值。

北航数值分析报告第一次大作业(幂法反幂法)

北航数值分析报告第一次大作业(幂法反幂法)

一、问题分析与算法描述1. 问题的提出:〔1〕用幂法、反幂法求矩阵的按摸最大和最小特征值,并求出相应的特征向量。

其中要求:迭代精度达到。

〔2〕用带双步位移的QR法求上述的全部特征值,并求出每一个实特征值相应的特征向量。

2. 算法的描述:(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。

其迭代格式为:终止迭代的控制选用。

幂法的使用条件为实矩阵A具有n个线性无关的特征向量,其相应的特征值满足不等式或幂法收敛速度与比值或有关,比值越小,收敛速度越快。

(2) 反幂法反幂法用于计算实矩阵A按摸最小的特征值,其迭代格式为:每迭代一次都要求解一次线性方程组。

当k足够大时,,可近似的作为矩阵A的属于的特征向量。

比值越小,收敛的越快。

反幂法要求矩阵A非奇异。

(3) 带双步位移的QR分解法QR方法适用于计算一般实矩阵的全部特征值,尤其适用于计算中小型实矩阵的全部特征值。

本算例中采用带双步位移的QR方法,可加速收敛,其迭代格式为:二、计算结果与分析1. 计算结果:(1) 幂法:初始条件:最大迭代次数L=1000;向量计算结果:第1次迭代结果:最大特征值:0.00000e+000第2次迭代结果:最大特征值:2.48910e+000 相对误差:1.00000e+000 第3次迭代结果:最大特征值:1.67719e+000 相对误差:第4次迭代结果:最大特征值:-2.10960e+000 相对误差:1.79503e+000 第5次迭代结果:最大特征值:-6.13203e-001 相对误差:2.44030e+000 ……第794次迭代结果:最大特征值:-1.97638e+000 相对误差:最大特征值:-1.97638e+000 相对误差:********************最终迭代结果***************特征值:-1.97638e+000 相对误差:迭代次数:795(2) 反幂法:初始条件:最大迭代次数L=1000;向量运行结果:第1次迭代结果:最大特征值:1.07542e+000第2次迭代结果:最大特征值:-3.66550e+000 相对误差:1.29339e+000 第3次迭代结果:最大特征值:1.22709e+001 相对误差:1.29871e+000 第4次迭代结果:最大特征值:-1.03421e+000 相对误差:1.28650e+001 第5次迭代结果:最大特征值:相对误差:……第995次迭代结果:最大特征值:相对误差:第996次迭代结果:最大特征值:相对误差:最大特征值:相对误差:第998次迭代结果:最大特征值:相对误差:第999次迭代结果:最大特征值:相对误差:第1000次迭代结果:最大特征值:相对误差:******************************超过最大设定迭代次数,迭代失败!(3) 带双步位移的QR法:初始条件:最大迭代次数L=1000;向量运行结果:全部特征值:特征向量〔经谱X数归一化〕:实特征值对应特征向量:-0.062705 -0.022368 0.304372 0.064466 0.521833 -0.157024 0.136942 -0.218108 0.250264 -0.043064 -0.228688 -0.184632 -0.072871 0.124721 0.029070 0.102566 -0.136358 0.167727 0.085747 0.546165 实特征值对应特征向量:-0.018001 0.019652 0.273447 0.070528 0.274896 -0.144015 0.048385 0.376439 -0.583051 -0.054008 -0.168682 -0.113430 -0.034709 0.009204 0.472291 0.125664 -0.190617 0.113145 0.046278 0.059871 实特征值对应特征向量:0.106861 0.087709 -0.024967 -0.020897 0.064302 0.034047 0.535143 0.046383 0.028832 0.003479-0.097276 -0.383801 0.089445 -0.039560 -0.036928 -0.021330 0.014811 0.705836 -0.108904 0.082022 实特征值对应特征向量:-0.055201 0.003399 0.242191 0.102847 0.372470 -0.372826 0.113953 0.240659 -0.310401 -0.076590 -0.244632 -0.192549 -0.077259 0.263328 0.201662 0.154166 -0.407814 0.186782 0.094649 0.173302 实特征值对应特征向量:0.427828 -0.546801 0.007822 -0.382580 0.025199 0.012788 0.033241 0.005389 -0.004065 0.043524 -0.032112 -0.044233 0.135395 -0.006564 0.001214 0.020165 0.011678 0.050001 -0.585765 0.013115 实特征值对应特征向量:0.236032 -0.139250 -0.008143 0.638527 -0.009049 -0.002911 -0.001307 0.003054 0.006515 -0.030134 0.012712 0.011368 -0.018792 -0.001753 -0.005749 -0.014290 -0.005292 -0.014591 0.717590 0.001369 实特征值对应特征向量:-0.227404 -0.048154 0.022615 0.297305 0.070372 0.039927 0.078503 0.015822 -0.012182 0.605334 -0.083616 -0.106270 -0.573963 -0.019907 0.003839 0.051362 0.036567 0.115613 0.332707 0.036954 实特征值对应特征向量:-0.027768 -0.051081 -0.159642 -0.054573 -0.084441 0.118378 0.029553 0.211088 0.203867 0.0486272. 结果分析以上三种方法中,幂法计算共进展了795次迭代才达到收敛,计算量较大,收敛性不好;反幂法计算结果未能收敛,通过进一步分析发现,这是因为反幂法迭代程序未考虑按模最小特征值为复数的情况,造成迭代失败。

北航研究生数值分析作业第二题

北航研究生数值分析作业第二题

北航研究生数值分析作业第二题北航研究生数值分析作业第二题:一、算法设计方案1.按照题目给出的矩阵定义对矩阵A赋初值:对应的函数为a_init();2.对矩阵A进行householder变换,使其拟上三角化:对应的函数为householder();3.输出拟上三角化后的A:对应的函数为aout(int);4.对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代(最大迭代次数L=500),逐个求出其特征值,对应的函数为eigen_a();中间包含两个子程序:calc_mk()和qr_analyze(),分别用来计算矩阵M k和对M k进行QR 分解并得到A k+1;5.输出QR分解过程完毕后的A及求得的特征向量:对应的函数为aout()和eigenvalout();6.对于在第三步中求得的每个实特征值,使用带原点平移的反幂法求出其对应的特征向量,对应的函数为eigenvec();其中包含一个解方程(A-μI)=y k-1的程序段。

这部分也用迭代完成,仍然将最大迭代次数L设置为500;7.输出矩阵A的特征向量,结束计算:对应的函数为eigenvecout()。

算法编译环境:vlsual c++6.0二、源程序如下:#include#include#define N 10 //矩阵阶数;#define EPSL 1.0e-12 //迭代的精度水平;#define L 500 //迭代最大次数;#define OUTPUTMODE 1 //输出格式:0--输出至屏幕,1--输出至文件double a[N][N], a2[N][N], eigen[N][N]; //声明矩阵A;double sa_re[N] = {0}, sa_im[N] = {0}; //声明矩阵的特征值数组;double u_init[N] = {2,1,2,1,2,1,2,1,2,1}; //定义反幂法中使用的初始向量u;//主程序开始;int main(){FILE *p;void a_init();void householder();void equal_zero(double matrix[N][N], int);void eigenvec();int eigen_a();void aout(int);void eigenvalout(int);void eigenvecout(int);if(OUTPUTMODE){p = fopen("Result.txt", "w+");fprintf(p, "计算结果:\n");fclose(p);}a_init(); //对矩阵A进行初始化;householder(); //对矩阵A进行拟上三角化;equal_zero(a, N); //对矩阵A的元素进行归零处理,消除误差;aout(OUTPUTMODE); //输出A;if(eigen_a()) printf("迭代超过最大次数,特征值求解结果可能不正确。

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

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

北京航空航天大学《数值分析》计算实习报告第一大题学 院:自动化科学与电气工程学院 专 业: 控制科学与工程 学 生 姓 名: 学 号: 教 师: 电 话: 完 成 日 期: 2015年11月6日北京航空航天大学Beijing University of Aeronautics and Astronautics实习题目:第一题 设有501501⨯的实对称矩阵A ,其中,064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。

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

2.求A 的与数4015011λλλμ-+=kk 最接近的特征值)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 ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤,可知绝对值最大特征值必为1λ和501λ其中之一,故可用幂法求出绝对值最大的特征值λ,如果λ=0,则1λ=λ,否则501λ=λ。

将矩阵A 进行一下平移:I -A A'λ= (1)对'A 用幂法求出其绝对值最大的特征值'λ,则A 的另一端点特征值1λ或501λ为'λ+λ。

s λ为按模最小特征值,||min ||5011i i s λλ≤≤=,可对A 使用反幂法求得。

北航数值分析-实习作业1(C语言详细注释)

北航数值分析-实习作业1(C语言详细注释)

《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a 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 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011,λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的(谱范数)条件数2)A (cond 和行列式A det要求1. 算法的设计方案(A 的所有零元素都不能存储)2. 全部源程序(详细注释)。

变量为double ,精度-1210=ε,输出为e 型12位有效数字3. 特征值s 5011,,λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值 4. 讨论迭代初始向量的选取对计算结果的影响,并说明原因解答:1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=,所以s λ是按模最小的特征值,直接运用反幂法可求得。

对于5011,λλ,一个是最大的特征值,一个是最小的特征值,不能确定两者的绝对值是否相等,因此必须首先假设||||5011λλ≠,然后运用幂法,看能否求得一个特征值,如果可以求得一个,证明A 是收敛的,求得的结果是正确的,然后对A 进行带原点平移的幂法,偏移量是前面求得的特征值,可以求得另一个特征值,最后比较这两个特征值,较大的特征值是501λ,较小的特征值就是1λ。

如果在假设的前提下,无法运用幂法求得按模最大的特征值,即此时A 不收敛,则需要将A 进行带原点平移的幂法,平移量可以选取1,再重复上述步骤即可求得两个特征值。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
相关文档
最新文档