约化对称矩阵为三对角对称矩阵
非对称特征值问题-基本概念
2.0003 0.8171 3.6516 , A25 0.0002 0.3336 3.7 263 0.7456 2.3333
2.0002 2.9999 2.2374 。 A26 0.0001 2.9996 2.23 66 2.2349 0.9998
我们称这种分块上三交阵为矩阵A的Schur分块上三角阵，
为了节省运算工作量，实用的方法是先将矩阵约化为与Schur分块上三角
阵相似的Hessenberg形。
（bij） Rnn的次对角线以下的元素bij =0 定义 若矩阵B
(i>j+1), 则称B为上Hessenberg矩阵，简称Hessenberg形，即
0.283205888 A3 0.157002612 0 0.287735078 A4 0.036401350 0
Pn2，则定理得证。
推论 对于任何对称矩阵A Rnn , 存在正交阵Q，使得B QT AQ为 对称三对角阵。
4.2
QR算法及其收敛性
QR算法可以用来求任意的非奇异矩阵的全部特征值，是目前计算这类问 题最有效的方法之一。它基于对任何实的非奇异矩阵都可以分解为正交阵Q和 上三角矩阵R的乘积。
设向量w Rn , w 2 1, 则称
H (w) I 2wwT
为（初等）镜面反射矩阵或Householder变换矩阵。
Houholder矩阵H=H（w）有如下性质：
(1) (2)
H是对称正交阵，即H H T H 1。
对任何x Rn , 记y Hx, 有 y 2 x 2
1 B
n 1
三种常用固有振动特征值解法的比较
2005全国结构动力学学术研讨会海南省海口市，2005.12.19-20中国振动工程学会结构动力学专业委员会三种常用固有振动特征值解法的比较宫玉才1周洪伟 陈 璞 袁明武（北京大学力学与工程科学系 北京，100871）Email ：yuanmw@摘要： 本文以高效的细胞稀疏直接快速解法为核心步骤，实现了快速的固有振动广义特征值问题解法, 并在相同的允许模态误差的意义下检验了三种结构动力学中常用的大型矩阵特征模态算法——子空间迭代法、迭代Ritz 向量法和迭代Lanczos 法的计算效率。
迭代Ritz 向量法平均而言最快，子空间迭代法最慢，三种解法效率相差不是太大。
与ANSYS 的子空间迭代和Lanczos 法相比，本文的子空间迭代比ANSYS 的效率高很多，Lanczos 法和ANSYS 的差不多 。
大量较大规模的例题显示，本文对特征值算法的改进是十分有效的，算法的健壮性，通用性都达到了高水平。
关键词：特征值，结构振动，迭代法，高效能计算1高等学校博士学科点专项科研基金资助项目 （编号：20030001112）引言在工程有限元分析中常常要求解广义代数特征值问题0K M ϕλϕ−= (1)的部分低阶特征值与特征向量。
对于矩阵阶数超过1000的大型问题，子空间迭代法、Ritz 向量法和Lanczos 法被公认为求解部分低阶极端特征值和特征向量的有效方法。
尽管国内外的有限元软件都提供广义代数特征值问题(1)的多种解法，但结果仍然不能令人完全满意，漏根与多根、自由模态误判都时有发生。
传统上，低端特征值问题求解过程极度依赖于谱变换的线性方程组()T K M x LDL x My µ−==(2)的解法，移轴矩阵K M µ−的LDLT 三角分解是计算量最大的主要步骤。
在以变带宽解法为核心步骤的特征值解法中，它常常占到特征值问题计算时间的70%到90%。
本文采用了文[1]提出的一个效率非常高的有限元解法-细胞稀疏直接快速解法（简称细胞解法）替换变带宽解法，极大地提高了三角分解的效率。
特征值
特征值与特征向量的计算与应用内容提要：很多实际问题可以归结为求矩阵的特征值问题，如结构动力学、电力网络、量子化学和物理学中某些临界值的确定等。
文章中给出了一些求解特征值和特征向量的计算方法。
包括幂法和反幂法，jocobi 和QR 算法和特征值与特征向量的应用。
关键词：特征值 特征向量 幂法 反幂法 QR 算法 特征值与特征向量的应用 一．引言为了利用矩阵研究线性变换，希望能找到线性空间的基使线性变换在该基下的矩阵具有最简单的形式，因此我们引进了特征值与特征向量。
特征值与特征向量在线性变换中起着举足轻重的作用，充分利用特征值与特征向量的命题与性质对我们解题带来极大的帮助，能使复杂的问题变得简单，化简为易，化繁为简。
定义：设A 是数域P 上线性空间V 的一个线性变换，如果对于数域P 中一数0λ，存在一个非零向量ξ，使得0A ξλξ=那么0λ称为A 的一个特征值，而ξ称为A 的属于特征值0λ的一个特征向量。
二．特征值与特征向量的计算方法1.幂法在实际问题中，常需要求矩阵的规模最大的特征值和相应的特征向量.对于这样的问题,幂法是合适的.幂法是计算主特征值的一种迭代方法,它的最大优点是方法简单,由于它主要是矩阵和向量的乘法运算,不改变原矩阵的稀疏结构,所以幂法对稀疏矩阵是适合的,但有时收敛速度很慢.设A 的主特征值为1λ,满足nλλλ≥≥> 21, (1)相应的n 个线性无关的特征向量为.,,,21n x x x任意非零的初始向量0v 可表示为这些特征向量的线性组合∑==ni i i x a v 10 ,其中01≠a ,构造迭代序列k k v A v=+1,k=0,1,…,则有(),111101k kn i i k i i k k k x a x a v A v A v ελλ+===∑=- (2)这里i kni i i k x a∑=⎪⎪⎭⎫⎝⎛=21λλε (3)从而()∞→→k x a v kk111λ.当11>λ时,k1λ增长非常快,计算k v 会导致溢出.为了解决这个问题,引入记号()i v v=max ,(4)其中i v 满足j nj i v v ≤≤=1max ,它表示取v 的按模最大的分量.构造迭代序列(称为幂法)()01110,0,,1,2,.max n i i i k k k k k v a x a u Av u v k u =-⎧=≠≠⎪⎪⎪=⎨⎪⎪==⎪⎩∑(5) 容易证明().2,1,0max =≠k u k (6)根据迭代上式有()().max max 11∏=-==ki ik k k k u v Ak u Av v（7）由于k v 的最大分量为1，从（7）式及（2）式得到()()110max max x x v A v A v kk k→=.另外,由于()101max v A v A Av u k k k k --==可知()()()()().max max max max max 1111111010λεελ→++==--k k k k k x a x a v A v A u设A 的主特征值为实重根,即121,,r r r n λλλλλλ+===>≥≥在上述推导中,把k ε修改为11kni i i i r a x λλ=+⎛⎫ ⎪⎝⎭∑ 及11a x 变为1ri i i a x =∑后，k v 还是会收敛到1λ对应的特征向量，即 ()()11111,max .max r rk k r r a x a x v u a x a x λ+→→+若特征值的分布为121,,r r r n λλλλλλ+===>≥≥ 则情况比较复杂，这里不做讨论.幂法的收敛速度和初始向量0v 的选取有关.当0v (不妨设021v =)在1x 上的投影1a 比较大,而在其他特征向量1x上的投影比较小时,根据(3)式知,此时幂法收敛会更快.这说明了当初是向量和主特征值对应的特征向量1x靠的很近时,幂法收敛会很快.幂法的收敛性主要由因子211λλ≈时,幂法收敛会很慢.此时可以用原点平移法进行加速,这可以看成使得收敛速度变快的一种预处理过程.设A 的特征值分布为12,n λλλ≥≥≥ (8)相应的特征向量为12,,n x x x.令I p A B-=,其中p 为待定系数，B的特征值分布为12,n p p p λλλ-≥-≥≥-B的特征向量和A 的特征向量相同.现取p 使得max ,j i i jp p λλ≠->-(9)则j p λ-为B的主特征值.对B施行幂法,求出它的主特征值j p λ-和特征向量j x ,从而也就知道了j λ和j x.对B 用幂法的收敛速度由因子max i i j jpp λλ≠--(10)确定.希望选取p 使得（10）式尽可能的小，我们把这种方法称为原点平移的幂法。
三对角矩阵的特征值
三对角矩阵的特征值一、引言矩阵在数学中有着广泛的应用，其中三对角矩阵是一类特殊的矩阵。
三对角矩阵是指除了主对角线和两条相邻的对角线之外，其余元素均为零的矩阵。
三对角矩阵在科学计算中有着广泛的应用，如求解常微分方程、偏微分方程、线性方程组等问题。
本文将详细介绍三对角矩阵的特征值及其求解方法。
二、三对角矩阵的特征值1. 定义设A是一个n×n实数或复数矩阵，则如果存在一个非零向量x使得Ax=λx，则称λ是A的特征值，x是A与λ相应的特征向量。
2. 三对角矩阵的特征值性质（1）三对角矩阵A只有n个特征值。
（2）若λ为A的一个特征值，则其代数重数等于几何重数。
（3）若λ为A的一个特征值，则它所属行列式|A-λE|中至少有一个元素不为零。
3. 三对角矩阵特征值求解方法（1）Jacobi法：Jacobi法是一种迭代求解特征值和特征向量的方法。
其基本思想是通过相似变换将矩阵对角化，使得对角线上的元素即为矩阵的特征值。
（2）QR方法：QR方法是一种迭代求解特征值和特征向量的方法。
其基本思想是通过正交相似变换将矩阵对称三对角化，然后通过反复施加Householder变换或Givens旋转来逐步将矩阵转化为上Hessenberg矩阵或上三角矩阵，最终得到矩阵的特征值。
（3）Lanczos算法：Lanczos算法是一种迭代求解特征值和特征向量的方法。
其基本思想是通过Krylov空间来逼近原始问题的最小或最大特征值，从而得到相应的特征向量。
三、Jacobi法求解三对角矩阵的特征值1. 算法流程（1）选取初始正交变换Q0=I。
（2）计算A1=QTAQ，其中Q为正交变换。
（3）寻找A1中非对角线元素中绝对值最大的元素aij，并确定旋转角度θ。
（4）构造Givens旋转矩阵G(i,j,θ)，使得G(i,j,θ)A1G(i,j,θ)T中元素a′ij=0。
（5）计算Q1=G(i,j,θ)Q0。
（6）将A2=Q1TA1Q1，重复步骤（3）至（5）直到A2为对角矩阵或满足一定的精度要求为止。
应用数值分析(第四版)课后习题答案第2章
第二章习题解答1.（1） R n×n中的子集“上三角阵”和“正交矩阵”对矩阵乘法是封闭的。
（2）R n×n中的子集“正交矩阵”，“非奇异的对称阵”和“单位上（下）三角阵”对矩阵求逆是封闭的。
设A 是ｎ×ｎ的正交矩阵。
证明A -1也是ｎ×ｎ的正交矩阵。
证明：(1),n nA B A B R⨯∈证明：为上三角阵，为上三角阵，10(),0(),0(),,()(()),()()ij ij nij ik kj ij k n n T T T T T T T T T T a i j b i j C AB c a b c i j A B A B R AA A A E BB B B EAB AB ABB A E AB AB B A AB E AB =⨯∴=>=>==∴=>∴∈========∴∑则上三角阵对矩阵乘法封闭。
以下证明：为正交矩阵，为正交矩阵，为正交矩阵,故正交矩阵对矩阵乘法封闭。
（2）A 是ｎ×ｎ的正交矩阵∴A A -1 =A -1A=E 故（A -1）-1=A∴A -1（A -1）-1=（A -1）-1A -1 =E 故A -1也是ｎ×ｎ的正交矩阵。
设A 是非奇异的对称阵，证A -1也是非奇异的对称阵。
A 非奇异 ∴A 可逆且A -1非奇异 又A T =A ∴（A -1）T =（A T ）-1=A-1故A -1也是非奇异的对称阵设A 是单位上（下）三角阵。
证A -1也是单位上（下）三角阵。
证明：A 是单位上三角阵，故|A|=1，∴A 可逆，即A -1存在，记为（b ij ）n×n由A A -1=E ，则∑==nj ik jkij ba 1δ （其中0=ij a j ＞i 时，1=ii a ）故b nn =1, b ni =0 (n≠j)类似可得，b ii =1 (j=1…n) b jk =0 (k ＞j)即A -1是单位上三角阵综上所述可得。
反幂法用来计算矩阵按模最小的特征值及其特征向量
a(k) kn
a(k) k 1,n
a(k) n,k 1
a(k) nn
17
其中
ck
(ak( k)1,k
,,
a(k) nk
)T
Rnk ,
A1(1k ) 为
k 阶上海森伯
格阵，A2(2k ) R . (nk )(nk )
设 ck 0，于是可选择初等反射阵 Rk 使 Rkck ke1， 其中, Rk 计算公式为
13
（1） 设
a11 a12 a1n
A
a21
an1
a22
an 2
a2n
ann
a11 c1
A(1) 12
A(1) 22
,
其中 c1 (a21,, an1)T R n1，不妨设 c1 0, 否则这一 步不需要约化.
本节讨论两个问题 （1） 用初等反射阵作正交相似变换约化一般实矩阵A 为上海森伯格阵. （2） 用初等反射阵作正交相似变换约化对称矩阵A 为 对称三对角阵. 于是，求原矩阵特征值问题，就转化为求上海森伯格阵 或对称三对角阵的特征值问题.
8.3.2 用正交相似变换约化一般矩阵为上海森柏格阵
设 A (aij ) R nn . 可选择初等反射阵U1,U2 ,,Un2 使 A 经正交相似变换约化为一个上海森伯格阵.
则
A2
U1 A1U1
a11 R1c1
A(1) 12
R1
R1
A(1) 22
R1
（3.1）
15
a11
数值分析第三课本习题及答案
第一章 绪 论1. 设x >0,x 的相对误差为δ,求ln x 的误差.2. 设x 的相对误差为2％,求nx 的相对误差.3. 以下各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指出它们是几位有效数字:*****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====⨯4. 利用公式(3.3)求以下各近似值的误差限:********12412324(),(),()/,i x x x ii x x x iii x x ++其中****1234,,,x x x x 均为第3题所给的数.5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少?6. 设028,Y =按递推公式1n n Y Y -=( n=1,2,…)计算到100Y .27.982(五位有效数字),试问计算100Y 将有多大误差?7. 求方程25610x x -+=的两个根,使它至少具有四位有效数字27.982).8. 当N 充分大时,怎样求211Ndx x +∞+⎰?9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2?10. 设212S gt =假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对误差增加,而相对误差却减小. 11. 序列{}n y 满足递推关系1101n n y y -=-(n=1,2,…),假设0 1.41y ≈(三位有效数字),计算到10y 时误差有多大?这个计算过程稳定吗?12.计算61)f =,1.4≈,利用以下等式计算,哪一个得到的结果最好?3--13.()ln(f x x =,求f (30)的值.假设开平方用六位函数表,问求对数时误差有多大?假设改用另一等价公式ln(ln(x x =-计算,求对数时误差有多大?14. 试用消元法解方程组{101012121010;2.x x x x +=+=假定只用三位数计算,问结果是否可靠?15. 已知三角形面积1sin ,2s ab c =其中c 为弧度,02c π<<,且测量a ,b ,c 的误差分别为,,.a b c ∆∆∆证明面积的误差s ∆满足.s a b cs a b c ∆∆∆∆≤++第二章 插值法1. 根据(2.2)定义的范德蒙行列式,令2000011211121()(,,,,)11n n n n n n n n n x x x V x V x x x x x x x xx x ----==证明()n V x 是n 次多项式,它的根是01,,n x x -,且101101()(,,,)()()n n n n V x V x x x x x x x ---=--.2. 当x = 1 , -1 , 2 时, f (x)= 0 , -3 , 4 ,求f (x )的二次插值多项式.3. 给出f (x )=ln x 的数值表用线性插值及二次插值计算ln 0.54 的近似值.4. 给出cos x ,0°≤x ≤90°的函数表,步长h =1′=(1/60)°,假设函数表具有5位有效数字,研究用线性插值求cos x 近似值时的总误差界.5. 设0k x x kh =+,k =0,1,2,3,求032max ()x x x l x ≤≤.6. 设jx 为互异节点(j =0,1,…,n ),求证:i)0()(0,1,,);nkkj jj x l x x k n =≡=∑ii)()()1,2,,).nk jj j xx l x k n =-≡0(=∑7. 设[]2(),f x C a b ∈且()()0f a f b ==,求证21()()().8max max a x ba xb f x b a f x ≤≤≤≤≤-"8. 在44x -≤≤上给出()xf x e =的等距节点函数表,假设用二次插值求xe 的近似值,要使截断误差不超过610-,问使用函数表的步长h 应取多少? 9. 假设2n n y =,求4n y ∆及4n y δ.10. 如果()f x 是m 次多项式,记()()()f x f x h f x ∆=+-,证明()f x 的k 阶差分()(0)k f x k m ∆≤≤是m k -次多项式,并且()0(m lf x l +∆=为正整数).11. 证明1()k k k k k k f g f g g f +∆=∆+∆.12. 证明110010.n n kkn n k k k k f gf g f g g f --+==∆=--∆∑∑13. 证明1200.n j n j y y y -=∆=∆-∆∑14. 假设1011()n n n n f x a a x a x a x --=++++有n 个不同实根12,,,n x x x ,证明{10,02;, 1.1()n k njk n a k n j jx f x -≤≤-=-=='∑15. 证明n 阶均差有以下性质: i)假设()()F x cf x =,则[][]0101,,,,,,n n F x x x cf x x x =;ii) 假设()()()F x f x g x =+,则[][][]010101,,,,,,,,,n n n F x x x f x x x g x x x =+.16. 74()31f x x x x =+++,求0172,2,,2f ⎡⎤⎣⎦及0182,2,,2f ⎡⎤⎣⎦.17. 证明两点三次埃尔米特插值余项是(4)22311()()()()/4!,(,)k k k k R x f x x x x x x ++=ξ--ξ∈并由此求出分段三次埃尔米特插值的误差限.18. 求一个次数不高于4次的多项式()P x ,使它满足(0)(1)P P k =-+并由此求出分段三次埃尔米特插值的误差限.19. 试求出一个最高次数不高于4次的函数多项式()P x ,以便使它能够满足以下边界条件(0)(0)0P P ='=,(1)(1)1P P ='=,(2)1P =.20. 设[](),f x C a b ∈,把[],a b 分为n 等分,试构造一个台阶形的零次分段插值函数()n x ϕ并证明当n →∞时,()n x ϕ在[],a b 上一致收敛到()f x .21. 设2()1/(1)f x x =+,在55x -≤≤上取10n =,按等距节点求分段线性插值函数()h I x ,计算各节点间中点处的()h I x 与()f x 的值,并估计误差.22. 求2()f x x =在[],a b 上的分段线性插值函数()h I x ,并估计误差.23. 求4()f x x =在[],a b 上的分段埃尔米特插值,并估计误差. 24. 给定数据表如下:试求三次样条插值并满足条件i) (0.25) 1.0000,(0.53)0.6868;S S '='=ii)(0.25)(0.53)0.S S "="=25. 假设[]2(),f x C a b ∈,()S x 是三次样条函数,证明 i)[][][][]222()()()()2()()()bbbba a a a f x dx S x dx f x S x dx S x f x S x dx "-"="-"+""-"⎰⎰⎰⎰;ii) 假设()()(0,1,,)i i f x S x i n ==,式中i x 为插值节点,且01n a x x x b =<<<=,则[][][]()()()()()()()()()baS x f x S x dx S b f b S b S a f a S a ""-"="'-'-"'-'⎰.26. 编出计算三次样条函数()S x 系数及其在插值节点中点的值的程序框图(()S x 可用(8.7)式的表达式).第三章 函数逼近与计算1. (a)利用区间变换推出区间为[],a b 的伯恩斯坦多项式.(b)对()sin f x x =在[]0,/2π上求1次和三次伯恩斯坦多项式并画出图形,并与相应的马克劳林级数部分和误差做比较. 2. 求证:(a)当()m f x M ≤≤时,(,)n m B f x M ≤≤. (b)当()f x x =时,(,)n B f x x =.3. 在次数不超过6的多项式中,求()sin 4f x x =在[]0,2π的最正确一致逼近多项式.4. 假设()f x 在[],a b 上连续,求()f x 的零次最正确一致逼近多项式.5. 选取常数a ,使301max x x ax≤≤-到达极小,又问这个解是否唯一?6. 求()sin f x x =在[]0,/2π上的最正确一次逼近多项式,并估计误差.7. 求()xf x e =在[]0,1上的最正确一次逼近多项式.8. 如何选取r ,使2()p x x r =+在[]1,1-上与零偏差最小?r 是否唯一? 9. 设43()31f x x x =+-,在[]0,1上求三次最正确逼近多项式. 10. 令[]()(21),0,1n n T x T x x =-∈,求***0123(),(),(),()T x T x T x T x .11. 试证{}*()nT x 是在[]0,1上带权ρ=的正交多项式.12. 在[]1,1-上利用插值极小化求11()f x tg x -=的三次近似最正确逼近多项式. 13. 设()xf x e =在[]1,1-上的插值极小化近似最正确逼近多项式为()n L x ,假设nf L ∞-有界,证明对任何1n ≥,存在常数n α、n β,使11()()()()(11).n n n n n T x f x L x T x x ++α≤-≤β-≤≤14. 设在[]1,1-上234511315165()128243843840x x x x x x ϕ=-----,试将()x ϕ降低到3次多项式并估计误差. 15. 在[]1,1-上利用幂级数项数求()sin f x x =的3次逼近多项式,使误差不超过0.005.16. ()f x 是[],a a -上的连续奇(偶)函数,证明不管n 是奇数或偶数,()f x 的最正确逼近多项式*()n n F x H ∈也是奇(偶)函数.17. 求a 、b 使[]220sin ax b x dxπ+-⎰为最小.并与1题及6题的一次逼近多项式误差作比较.18. ()f x 、[]1(),g x C a b ∈,定义 ()(,)()();()(,)()()()();b baaa f g f x g x dxb f g f x g x dx f a g a =''=''+⎰⎰问它们是否构成内积?19. 用许瓦兹不等式(4.5)估计6101x dx x +⎰的上界,并用积分中值定理估计同一积分的上下界,并比较其结果.20. 选择a ,使以下积分取得最小值:1122211(),x ax dx x ax dx----⎰⎰.21. 设空间{}{}10010121,,,span x span x x 1ϕ=ϕ=,分别在1ϕ、2ϕ上求出一个元素,使得其为[]20,1x C ∈的最正确平方逼近,并比较其结果.22. ()f x x =在[]1,1-上,求在{}2411,,span x x ϕ=上的最正确平方逼近.23.sin (1)arccos ()n n x u x +=是第二类切比雪夫多项式,证明它有递推关系()()()112n n n u x xu x u x +-=-.24. 将1()sin 2f x x=在[]1,1-上按勒让德多项式及切比雪夫多项式展开,求三次最正确平方逼近多项式并画出误差图形,再计算均方误差.25. 把()arccos f x x =在[]1,1-上展成切比雪夫级数.26. 用最小二乘法求一个形如2y a bx =+的经验公式,使它与以下数据拟合,并求均方误差.27.28. 在某化学反应里,根据实验所得分解物的浓度与时间关系如下:用最小二乘拟合求.29. 编出用正交多项式做最小二乘拟合的程序框图. 30. 编出改良FFT 算法的程序框图. 31. 现给出一张记录{}{}4,3,2,1,0,1,2,3k x =,试用改良FFT 算法求出序列{}k x 的离散频谱{}k C (0,1,,7).k =第四章 数值积分与数值微分1. 确定以下求积公式中的待定参数,使其代数精度尽量高,并指明所构造出的求积公式所具有的代数精度: (1)101()()(0)()hh f x dx A f h A f A f h --≈-++⎰; (2)21012()()(0)()hh f x dx A f h A f A f h--≈-++⎰;(3)[]1121()(1)2()3()/3f x dx f f xf x -≈-++⎰;(4)[][]20()(0)()/1(0)()hf x dx h f f h ah f f h ≈++'-'⎰.2. 分别用梯形公式和辛普森公式计算以下积分:(1)120,84xdx n x =+⎰; (2)1210(1),10x e dx n x --=⎰;(3)1,4n =⎰; (4),6n =.3. 直接验证柯特斯公式(2.4)具有5次代数精度.4. 用辛普森公式求积分1x e dx-⎰并计算误差.5. 推导以下三种矩形求积公式:(1)2()()()()()2ba f f x dxb a f a b a 'η=-+-⎰; (2)2()()()()()2baf f x dx b a f b b a 'η=---⎰;(3)3()()()()()224baa b f f x dx b a f b a +"η=-+-⎰.6. 证明梯形公式(2.9)和辛普森公式(2.11)当n →∞时收敛到积分()baf x dx⎰.7.用复化梯形公式求积分()baf x dx⎰,问要将积分区间[],a b 分成多少等分,才能保证误差不超过ε(设不计舍入误差)?8.1x e dx-,要求误差不超过510-.9. 卫星轨道是一个椭圆,椭圆周长的计算公式是S a =θ,这里a 是椭圆的半长轴,c是地球中心与轨道中心(椭圆中心)的距离,记h 为近地点距离,H 为远地点距离,6371R =公里为地球半径,则(2)/2,()/2a R H h c H h =++=-.我国第一颗人造卫星近地点距离439h =公里,远地点距离2384H =公里,试求卫星轨道的周长.10. 证明等式3524sin3!5!n nn n ππππ=-+-试依据sin(/)(3,6,12)n n n π=的值,用外推算法求π的近似值.11. 用以下方法计算积分31dyy ⎰并比较结果.(1) 龙贝格方法;(2) 三点及五点高斯公式;(3) 将积分区间分为四等分,用复化两点高斯公式.12. 用三点公式和五点公式分别求21()(1)f x x =+在x =1.0,1.1和1.2处的导数值,并估计误差.()f x 的值由下表给出:第五章 常微分方程数值解法1. 就初值问题0)0(,=+='y b ax y 分别导出尤拉方法和改良的尤拉方法的近似解的表达式，并与准确解bx ax y +=221相比较。
Lanczos算法
Lanczos算法是一种将对称矩阵通过正交相似变换变成对称三对角矩阵的算法，以20世纪匈牙利数学家Cornelius Lanczos命名。
Lanczos算法实际上是Arnoldi算法对于对称矩阵的特殊形式，可应用于对称矩阵线性方程组求解的Krylov子空间方法以及对称矩阵的特征值问题。
Lanczos算法给定对称矩阵A;选取单位向量v_1;设定v_0为零向量;设定b_0=0;for i=1:ma_i=(Av_i,v_i);b_i=||Av_i-a_iv_i-b_{i-1}v_{i-1}||;b_i v_{i+1} = Av_i - a_i v_i - b_{i-1}v_{i-1};end由上述Lanczos算法得：V'AV=T,其中V=[v_1,...,v_m], T=tridiag(b,a,b), a=[a_1,...,a_m],b=[b_1,...,b_m].编辑本段该算法的一个matlab实现程序A代表任意一个需要三对角化的矩阵，b是任意一个向量，且b的行数与A的列数相同因为要用到v = A*q;nmax是你想要得到的矩阵的大小，例如nmax=12，最后得到12*12的三对角矩阵。
结果输出的是一个三对角矩阵输入形式为：lanczos([1 2 3;4 5 6;7 8 9],[1;1;1],12);function T = lanczos(A, b, nmax)m = size(A,1);beta(1) = 0;qprev = zeros(m, 1);q = b / norm(b);for n = 1:nmaxv = A*q;alpha(n) = q' * v;v = v - beta(n) * qprev - alpha(n) * q;beta(n+1) = norm(v);qprev = q;q = v / beta(n+1);endbeta = beta(2:end-1);T = diag(alpha) + diag(beta,1) + diag(beta,-1);。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
//约化对称矩阵为三对角对称矩阵//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵//a-长度为n*n的数组，存放n阶实对称矩阵//n-矩阵的阶数//q-长度为n*n的数组，返回时存放Householder变换矩阵//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素void eastrq(double a[],int n,double q[],double b[],double c[]);////////////////////////////////////////////////////////////////求实对称三对角对称矩阵的全部特征值及特征向量//利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量//n-矩阵的阶数//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素//q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组// 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组//a-长度为n*n的数组，存放n阶实对称矩阵int ebstq(int n,double b[],double c[],double q[],double eps,int l);////////////////////////////////////////////////////////////////约化实矩阵为赫申伯格(Hessen berg)矩阵//利用初等相似变换将n阶实矩阵约化为上H矩阵//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵//n-矩阵的阶数void echbg(double a[],int n);////////////////////////////////////////////////////////////////求赫申伯格(Hessen berg)矩阵的全部特征值//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//利用带原点位移的双重步QR方法求上H矩阵的全部特征值//a-长度为n*n的数组，存放上H矩阵//n-矩阵的阶数//u-长度为n的数组，返回n个特征值的实部//v-长度为n的数组，返回n个特征值的虚部//eps-控制精度要求//jt-整型变量，控制最大迭代次数int edqr(double a[],int n,double u[],double v[],double eps,int jt);////////////////////////////////////////////////////////////////求实对称矩阵的特征值及特征向量的雅格比法//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值//n-矩阵的阶数//u-长度为n*n的数组，返回特征向量(按列存储)//eps-控制精度要求//jt-整型变量，控制最大迭代次数int eejcb(double a[],int n,double v[],double eps,int jt);//////////////////////////////////////////////////////////////选自<<徐世良数值计算程序集(C)>>每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。
今天都给贴出来了#include "stdio.h"#include "math.h"//约化对称矩阵为三对角对称矩阵//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵//a-长度为n*n的数组，存放n阶实对称矩阵//n-矩阵的阶数//q-长度为n*n的数组，返回时存放Householder变换矩阵//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素void eastrq(double a[],int n,double q[],double b[],double c[]){int i,j,k,u,v;double h,f,g,h2;for (i=0; i<=n-1; i++){for (j=0; j<=n-1; j++){u=i*n+j; q[u]=a[u];}}for (i=n-1; i>=1; i--){h=0.0;if (i>1){for (k=0; k<=i-1; k++){u=i*n+k;h=h+q[u]*q[u];}}if (h+1.0==1.0){c[i-1]=0.0;if (i==1){c[i-1]=q[i*n+i-1];}b[i]=0.0;}else{c[i-1]=sqrt(h);u=i*n+i-1;if (q[u]>0.0){c[i-1]=-c[i-1];}h=h-q[u]*c[i-1];q[u]=q[u]-c[i-1];f=0.0;for (j=0; j<=i-1; j++) {q[j*n+i]=q[i*n+j]/h;g=0.0;for (k=0; k<=j; k++) {g=g+q[j*n+k]*q[i*n+k];}if (j+1<=i-1){for (k=j+1; k<=i-1; k++){g=g+q[k*n+j]*q[i*n+k];}}c[j-1]=g/h;f=f+g*q[j*n+i];}h2=f/(h+h);for (j=0; j<=i-1; j++) {f=q[i*n+j];g=c[j-1]-h2*f;c[j-1]=g;for (k=0; k<=j; k++){u=j*n+k;q[u]=q[u]-f*c[k-1]-g*q[i*n+k];}}b[i]=h;}}b[0]=0.0;for (i=0; i<=n-1; i++){if ((b[i]!=0.0)&&(i-1>=0)){for (j=0; j<=i-1; j++){g=0.0;for (k=0; k<=i-1; k++){g=g+q[i*n+k]*q[k*n+j];}for (k=0; k<=i-1; k++){u=k*n+j;q[u]=q[u]-g*q[k*n+i];}}}u=i*n+i;b[i]=q[u];q[u]=1.0;if (i-1>=0){for (j=0; j<=i-1; j++){q[i*n+j]=0.0;q[j*n+i]=0.0;}}}return;//求实对称三对角对称矩阵的全部特征值及特征向量//利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量//n-矩阵的阶数//b-长度为n的数组，返回时存放三对角阵的主对角线元素//c-长度为n的数组，返回时前n-1个元素存放次对角线元素//q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组// 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组//a-长度为n*n的数组，存放n阶实对称矩阵int ebstq(int n,double b[],double c[],double q[],double eps,int l){int i,j,k,m,it,u,v;double d,f,h,g,p,r,e,s;c[n-1]=0.0;d=0.0;f=0.0;for (j=0; j<=n-1; j++){it=0;h=eps*(fabs(b[j])+fabs(c[j]));if (h>d){d=h;}m=j;while ((m<=n-1)&&(fabs(c[m])>d)){m=m+1;}if (m!=j){do{if (it==l){printf("fail\n");return(-1);}it=it+1;g=b[j];p=(b[j+1]-g)/(2.0*c[j]);if (p>=0.0){b[j]=c[j]/(p+r);}else{b[j]=c[j]/(p-r);}h=g-b[j];for (i=j+1; i<=n-1; i++){b[i]=b[i]-h;}f=f+h;p=b[m];e=1.0;s=0.0;for (i=m-1; i>=j; i--){g=e*c[i];h=e*p;if (fabs(p)>=fabs(c[i])) {e=c[i]/p;r=sqrt(e*e+1.0);c[i+1]=s*p*r;s=e/r;e=1.0/r;}else{e=p/c[i];r=sqrt(e*e+1.0);c[i+1]=s*c[i]*r;s=1.0/r;e=e/r;}p=e*b[i]-s*g;b[i+1]=h+s*(e*g+s*b[i]);for (k=0; k<=n-1; k++) {u=k*n+i+1;v=u-1;q[u]=s*q[v]+e*h;q[v]=e*q[v]-s*h;}}c[j]=s*p;b[j]=e*p;}while (fabs(c[j])>d);}b[j]=b[j]+f;}for (i=0; i<=n-1; i++){k=i; p=b[i];if (i+1<=n-1){j=i+1;while ((j<=n-1)&&(b[j]<=p)) {k=j;p=b[j];j=j+1;}}if (k!=i){b[k]=b[i];b[i]=p;for (j=0; j<=n-1; j++){u=j*n+i;v=j*n+k;p=q[u];q[u]=q[v];q[v]=p;}}}return(1);}//约化实矩阵为赫申伯格(Hessen berg)矩阵//利用初等相似变换将n阶实矩阵约化为上H矩阵//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵//n-矩阵的阶数void echbg(double a[],int n){ int i,j,k,u,v;double d,t;for (k=1; k<=n-2; k++){d=0.0;for (j=k; j<=n-1; j++){u=j*n+k-1;t=a[u];if (fabs(t)>fabs(d)){d=t;i=j;}}if (fabs(d)+1.0!=1.0){if (i!=k){for (j=k-1; j<=n-1; j++){u=i*n+j;v=k*n+j;t=a[u];a[u]=a[v];a[v]=t;}for (j=0; j<=n-1; j++){u=j*n+i;v=j*n+k;t=a[u];a[u]=a[v];a[v]=t;}}for (i=k+1; i<=n-1; i++){u=i*n+k-1;t=a[u]/d;a[u]=0.0;for (j=k; j<=n-1; j++){v=i*n+j;a[v]=a[v]-t*a[k*n+j];}for (j=0; j<=n-1; j++){v=j*n+k;a[v]=a[v]+t*a[j*n+i];}}}}return;}//求赫申伯格(Hessen berg)矩阵的全部特征值//利用带原点位移的双重步QR方法求上H矩阵的全部特征值//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//a-长度为n*n的数组，存放上H矩阵//n-矩阵的阶数//u-长度为n的数组，返回n个特征值的实部//v-长度为n的数组，返回n个特征值的虚部//eps-控制精度要求//jt-整型变量，控制最大迭代次数int edqr(double a[],int n,double u[],double v[],double eps,int jt) {int m,it,i,j,k,l,ii,jj,kk,ll;double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;it=0;m=n;while (m!=0){l=m-1;while((l>0)&&(fabs(a[l*n+l-1])>eps*(fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l])))) {l=l-1;}ii=(m-1)*n+m-1;jj=(m-1)*n+m-2;kk=(m-2)*n+m-1;ll=(m-2)*n+m-2;if (l==m-1){u[m-1]=a[(m-1)*n+m-1];v[m-1]=0.0;m=m-1; it=0;}else if (l==m-2){b=-(a[ii]+a[ll]);c=a[ii]*a[ll]-a[jj]*a[kk];w=b*b-4.0*c;y=sqrt(fabs(w));if (w>0.0){xy=1.0;if (b<0.0){xy=-1.0;}u[m-1]=(-b-xy*y)/2.0;u[m-2]=c/u[m-1];v[m-1]=0.0; v[m-2]=0.0; }else{u[m-1]=-b/2.0;u[m-2]=u[m-1];v[m-1]=y/2.0;v[m-2]=-v[m-1];}m=m-2;it=0;}else{if (it>=jt){printf("fail\n");return(-1);}it=it+1;for (j=l+2; j<=m-1; j++){a[j*n+j-2]=0.0;}for (j=l+3; j<=m-1; j++){a[j*n+j-3]=0.0;}for (k=l; k<=m-2; k++){if (k!=l){p=a[k*n+k-1];q=a[(k+1)*n+k-1];r=0.0;if (k!=m-2){r=a[(k+2)*n+k-1];}}else{x=a[ii]+a[ll];y=a[ll]*a[ii]-a[kk]*a[jj];ii=l*n+l;jj=l*n+l+1;kk=(l+1)*n+l;ll=(l+1)*n+l+1;p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y;q=a[kk]*(a[ii]+a[ll]-x);r=a[kk]*a[(l+2)*n+l+1];}if ((fabs(p)+fabs(q)+fabs(r))!=0.0) {xy=1.0;if (p<0.0){xy=-1.0;}s=xy*sqrt(p*p+q*q+r*r);if (k!=l){a[k*n+k-1]=-s;}e=-q/s;x=-p/s;y=-x-f*r/(p+s);g=e*r/(p+s);z=-x-e*q/(p+s);for (j=k; j<=m-1; j++) {ii=k*n+j;jj=(k+1)*n+j;p=x*a[ii]+e*a[jj];q=e*a[ii]+y*a[jj];r=f*a[ii]+g*a[jj];if (k!=m-2){kk=(k+2)*n+j;p=p+f*a[kk];q=q+g*a[kk];r=r+z*a[kk]; a[kk]=r;}a[jj]=q;a[ii]=p;}j=k+3;if (j>=m-1){j=m-1;}for (i=l; i<=j; i++){ii=i*n+k;jj=i*n+k+1;p=x*a[ii]+e*a[jj];q=e*a[ii]+y*a[jj];r=f*a[ii]+g*a[jj];if (k!=m-2){kk=i*n+k+2;p=p+f*a[kk];q=q+g*a[kk];r=r+z*a[kk]; a[kk]=r;}a[jj]=q;}}}}}return(1);}//求实对称矩阵的特征值及特征向量的雅格比法//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量//返回值小于0表示超过迭代jt次仍未达到精度要求//返回值大于0表示正常返回//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值//n-矩阵的阶数//u-长度为n*n的数组，返回特征向量(按列存储)//eps-控制精度要求//jt-整型变量，控制最大迭代次数int eejcb(double a[],int n,double v[],double eps,int jt){int i,j,p,q,u,w,t,s,l;double fm,cn,sn,omega,x,y,d;l=1;for (i=0; i<=n-1; i++){v[i*n+i]=1.0;for (j=0; j<=n-1; j++){if (i!=j){v[i*n+j]=0.0;}}}while (1==1){fm=0.0;for (i=0; i<=n-1; i++){for (j=0; j<=n-1; j++){d=fabs(a[i*n+j]);if ((i!=j)&&(d>fm)){fm=d;p=i;q=j;}}}if (fm<eps){return(1);}if (l>jt){return(-1);}l=l+1;u=p*n+q;w=p*n+p;t=q*n+p;s=q*n+q;x=-a[u];y=(a[s]-a[w])/2.0;omega=x/sqrt(x*x+y*y);if (y<0.0){omega=-omega;}sn=1.0+sqrt(1.0-omega*omega);sn=omega/sqrt(2.0*sn);cn=sqrt(1.0-sn*sn);fm=a[w];a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;a[u]=0.0;a[t]=0.0;for (j=0; j<=n-1; j++){if ((j!=p)&&(j!=q)){u=p*n+j;w=q*n+j;fm=a[u];a[u]=fm*cn+a[w]*sn;a[w]=-fm*sn+a[w]*cn;}}for (i=0; i<=n-1; i++){if ((i!=p)&&(i!=q)){u=i*n+p;w=i*n+q;fm=a[u];a[u]=fm*cn+a[w]*sn;a[w]=-fm*sn+a[w]*cn;}}for (i=0; i<=n-1; i++){u=i*n+p;w=i*n+q;fm=v[u];v[u]=fm*cn+v[w]*sn;v[w]=-fm*sn+v[w]*cn; }}return(1);}。