数值分析报告.doc
《数值分析报告.doc》由会员分享,可在线阅读,更多相关《数值分析报告.doc(21页珍藏版)》请在冰点文库上搜索。
北航2008级研究生《数值分析A》计算实习题
2008年11月15日
目录
1引言 1
2算法设计方案 1
3特殊情况处理 4
4计算结果 4
5结论 7
6参考文献 7
7附录 8
1引言
算法背景:
幂法,Jacobi方法及QR方法是求矩阵特征值和特征向量的常用数值方法,它们都是造构造迭代产生的矩阵序列来达到目的的。
幂法计算简单,特别适用于高阶稀疏矩阵,但其收敛速度不能令人满意,要想加快幂法的收敛速度可采用反幂法及位移技术。
Jacobi方法是古典方法,它收敛快,精度高,便于并行计算且算法稳定。
用Jacobi方法求出的特征向量在较好的正交性,不过它的计算量较大,当阶数n增大时收敛速度减慢,因此Jacobi方法适用于求低阶的对称矩阵的全部特征值和特征向量。
由J.G.Francis提出的QR算法的基本思想源于LR算法,即对矩阵分解获得一个相似于原矩阵的序列,使其收敛到一个易于求得特征值的形式。
LR算法比QR算法收敛速度快,但不稳定。
QR方法是60年代发展起来的,被人们称为数值数学,最值得注意的算法之一,它是目前求任意矩阵全部特征值和特征向量最有效的方法。
利用矩阵的QR分解,通过逆序相乘产生对原矩阵的一系列正交相似变换,使其变化为一个近似的上三角矩阵来求全部特征值。
这里QR分解是指将矩阵化为一个正交矩阵Q和一个上三角矩阵左乘的形式。
但是基本QR算法的收敛往往过慢,因此通常采用带原点位移的QR算法。
又由于A是实矩阵,它可能有复特征值,如果采用一般带原点位移的QR算法,变换中带有复位移量,则造成复数运算。
故而采用带双步位移的QR算法,可以减少迭代次数加速收敛,避免复数运算,在计算上是比较经济的。
2算法设计方案
2.1题目
试用带双步位移的QR分解法求矩阵A=的全部特征值,并对其中的每一个实特征值求相应的特征向量。
要求程序输出矩阵A经过拟上三角化后所得到的矩阵;对矩阵进行QR分解方法结束后所得的矩阵;矩阵A的全部特征值及相应于实特征值的特征向量。
采用e型输出数据,并且至少显示12位有效数字。
要求迭代的精度水平为=。
已知
2.2算法
2.2.1用带双步位移的QR方法求实矩阵全部特征值的算法
(1)使用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵;给定精度水平和迭代最大次数L。
【矩阵A的拟上三角化的算法如下:
记A,并记。
对于r=1,2,,n-2执行
(a)若全为零,则令,转(e);否则转(b)。
(b)计算,(若,则取),
(c)令
(d)计算,,,,
(e)继续】
(2)记,令k=1,m=n。
(3)如果,则得到A的一个特征值,置m=n-1,转(4);否则转(5)。
(4)如果m=1,则得到A的一个特征值,转(11);如果m=0,则直接转(11);如果m1,则转(3)。
(5)求二阶子阵的两个特征值和,即计算二次方程的两个根和。
(6)如果m=2,则得到A的两个特征值和,转(11);否则转(7)。
(7)如果,则得到A的两个特征值和,置m=m-2,转(4);否则转(8)。
(8)如果k=L,则计算终止,未得到A的全部特征值,否则转(9)。
(9)记,计算
(为m阶单位矩阵)
(对作QR分解)
【对作QR分解与的计算算法如下:
记,,。
对于r=1,2,,m-1执行
(a)若全为零,则令,转(e);否则转(b)。
(b)计算,(若,则取),
(c)令
(d)计算,,,,,,
(e)继续】
(10)置k=k+1,转(3)。
(11)A的全部特征值以计算完毕,停止计算。
2.2.2列主元素Gauss消去法求特征值对应的特征向量的算法
列主元素Gauss消去法具体算法如下:
记,
1.消元过程
对于1,2,,n-1执行
(1)选行号,使
(2)交换与以及与所含的数值
(3)对于i=k+1,k+2,,n计算
,
2.回代过程
3特殊情况处理
(1)由于拟上三角化和QR分解的算法有相同的部分故将两部分的子程序集成在一起(当
m=1时为矩阵的拟上三角化程序和当m=0时为QR分解程序)
(2)当输出“未求出矩阵A的全部特征值”时,检查L的大小,增大L的数值再计算。
4计算结果
4.1矩阵A经过拟上三角化后得到的矩阵
(由于一行放不下十个元素,故将矩阵每一行的十个元素分成两行列出)
矩阵A经过拟上三角化后所得的矩阵A(n-1):
-8.82751675883e-001 -9.93313649183e-002 -1.10334928599e+000 -7.60044358564e-001 1.54910107991e-001
-1.94659186287e+000 -8.78243638293e-002 -9.25588938718e-001 6.03259944053e-001 1.51886095647e-001
-2.34787836242e+000 2.37237010494e+000 1.81929082221e+000 3.23780410155e-001 2.20579844032e-001
2.10269266255e+000 1.81613808610e-001 1.27883908999e+000 -6.38057812440e-001 -4.15407560380e-001
-1.05568548241e-016 1.72827459997e+000 -1.17146764279e+000 -1.24383926270e+000 -6.39975834174e-001
-2.00283307904e+000 2.92494720612e-001 -6.41283006839e-001 9.78399762128e-002 2.55776357416e-001
-5.39338381277e-017 -3.16587386518e-017 -1.29166953413e+000 -1.11160351340e+000 1.17134682410e+000
-1.30735603002e+000 1.80369917775e-001 -4.24638535837e-001 7.98895523930e-002 1.60881992807e-001
1.53346499662e-017 5.96332140618e-017 -3.27306215536e-017 1.56012629853e+000 8.12504939751e-001
4.42175683292e-001 -3.58861612814e-002 4.69174231367e-001 -2.73659505009e-001 -7.35933465775e-002
1.30056273723e-016 -3.09706001089e-017 8.55612756587e-017 0.00000000000e+000 -7.70777375519e-001
-1.58305142574e+000 -3.04284317680e-001 2.52871244603e-001 -6.70992540145e-001 2.54461992908e-001
1.61021672477e-016 -2.21157183737e-016 -3.92548490401e-017 0.00000000000e+000 0.00000000000e+000
-7.46345345694e-001 -2.70836515702e-002 -9.48652189368e-001 1.19587108150e-001 1.92926561795e-002
1.36855018620e-016 7.15151319080e-017 -8.67261531120e-017 0.00000000000e+000 0.00000000000e+000
-1.07207471836e-016 -7.70180137436e-001 -4.69762399062e-001 4.98825946801e-001 1.13769160378e-001
-2.78085130072e-017 -6.70863078836e-017 8.47612717394e-017 0.00000000000e+000 0.00000000000e+000
-4.87977428748e-017 1.85482928622e-016 7.01316709211e-001 1.58218068848e-001 3.86259461423e-001
-2.12460444005e-017 -1.70797975893e-016 -2.79894228707e-017 0.00000000000e+000 0.00000000000e+000
4.15375832524e-017 1.22262169189e-016 0.00000000000e+000 4.84380760278e-001 3.99277799518e-001
4.2对矩阵进行QR分解后所得的矩阵
(由于一行放不下十个元素,故将矩阵每一行的十个元素分成两行列出)
对拟三角化后的矩阵A(n-1)进行QR分解方法结束后的矩阵:
1.53147700891e+000 -2.31151347770e+000 -5.13352162117e-002 -1.56228767543e-001 3.95270218975e-001
-4.60307123730e-001 -4.92665761043e-001 1.54615717231e-001 4.58107069729e-002 9.01038975186e-002
-3.19420475749e+000 -5.63497701019e-001 -2.71483239643e-001 -6.42858659429e-001 1.03913640971e+000
2.87072983195e+000 -7.00811288884e-001 -8.03485266055e-001 -7.59721306153e-001 2.13060816678e+000
3.55229096910e-016 1.78505809730e+000 -1.74399102792e+000 7.97374725643e-001 -6.25624360652e-001
-1.11600242269e+000 2.59559452545e-001 3.42921412680e-001 5.10113644468e-001 -9.75562845054e-001
-2.38895193340e-016 -2.87790678694e-016 -4.06561715572e-001 -5.06868432472e-001 1.82855102452e+000
5.18941002166e-001 -3.59510952371e-001 -1.13845216101e-001 -4.51225605836e-002 5.59740463210e-001
1.13201004533e-016 -9.76111965681e-018 3.31049963415e-016 1.48676607216e+000 3.71340891995e-001
-1.69622855201e-001 1.74920002091e-001 -1.73910145562e-001 -2.43180817397e-001 1.67045385533e-001
-1.95011078163e-017 -7.21884669449e-017 -8.32462282641e-017 -4.58793017954e-017 -6.55122789790e-001
-9.39665990713e-001 2.83922343339e-001 -5.65703210957e-002 -8.96845791791e-002 -3.88998224402e-001
1.55950784591e-016 1.55024068510e-016 1.39811287854e-016 -4.23319394901e-017 -2.36039665905e-016
-4.57356256434e-001 -1.25312729615e+000 -2.30116600863e-001 2.99531681862e-001 -4.26239601973e-001
4.17370606398e-017 4.99418388088e-017 4.85013749449e-017 -1.25802674480e-017 -1.44155149365e-016
-6.10890682971e-017 -4.66929122722e-001 8.08204282850e-001 -1.55343840859e-001 1.63069004449e-001
2.43149869404e-017 4.06334732673e-017 2.38568246195e-017 -1.93516231771e-017 -7.96194310258e-017
-3.33232735828e-017 -1.35991797860e-016 -1.59398041636e-001 9.12617646163e-002 3.11904439743e-003
-2.01582627298e-017 -5.01986969674e-017 -1.50289426382e-018 1.07390932683e-017 -8.37186191687e-017
-3.67783217967e-017 2.76682192903e-017 -5.82867087928e-016 -2.41769845651e-001 7.01517104511e-001
4.5矩阵A的全部特征值
矩阵A的全部特征值:
6.36062787575e-001 +0.00000000000e+000i
5.65048899350e-002 +0.00000000000e+000i
9.35588907819e-001 +0.00000000000e+000i
-9.80530956290e-001 +1.13948912743e-001i
-9.80530956290e-001 -1.13948912743e-001i
-1.48403982226e+000 +0.00000000000e+000i
1.57754855711e+000 +0.00000000000e+000i
-2.32349621021e+000 +8.93040517720e-001i
-2.32349621021e+000 -8.93040517720e-001i
3.38303961744e+000 +0.00000000000e+000i
4.6矩阵A的相应于实特征值的特征向量
(由于一行放不下十个元素,故将矩阵每一行的十个元素分成两行列出)
特征值对应的特征向量(按行排列)(不是实特征值的,不求对应的特征向量):
-1.07037468397e-001 -7.12345583256e-002 -3.90237917138e-001 4.46655513756e-002 7.19034744275e-001
-1.75815688440e-001 2.26537961544e-001 -3.76886150502e-001 -2.95625408505e-001 -2.25577972559e-002
2.09017518582e-001 2.00063796685e-001 -3.89176061817e-001 2.77939118707e-002 3.93236547906e-001
1.24703810952e-001 -6.44810939558e-001 3.02779853108e-001 2.91095263314e-001 -4.09436555812e-002
-8.05651548269e-002 -4.61113468803e-002 1.50243791976e-002 4.81208346184e-002 3.53633892069e-001
-2.08919077811e-001 1.55745966395e-001 -8.19670427271e-001 3.50982512819e-001 -2.88513852783e-002
不是实特征值,不求对应的特征向量
不是实特征值,不求对应的特征向量
5.60118116800e-001 -7.79341508770e-001 -1.33780114857e-002 2.77409287891e-001 -3.00557588301e-003
2.53483408960e-003 2.06284849088e-002 1.10134829107e-002 1.22486166525e-002 -3.23620930409e-002
-6.24962559917e-002 1.12090045150e-002 2.49666705179e-001 1.31364486150e-001 3.83541256120e-001
-8.15944310353e-001 1.24560352067e-001 6.83561054721e-002 -2.70587401777e-001 -1.00519108146e-001
不是实特征值,不求对应的特征向量
不是实特征值,不求对应的特征向量
1.05221573436e-001 2.18364197327e-001 4.73017877776e-001 2.60887478870e-001 3.05805625140e-001
2.58379783221e-001 -8.73379363951e-002 -4.05454033853e-001 -5.09013113720e-001 -2.40925044562e-001
5结论
QR方法是目前求任意矩阵全部特征值和特征向量最有效的方法。
利用矩阵的QR分解,通过逆序相乘产生对原矩阵的一系列正交相似变换,使其变化为一个近似的上三角矩阵(拟上三角矩阵)来求全部特征值,算出特征值后对实特征值利用列主元素Gauss消去法求实特征值对应的特征向量。
利用QR分解将矩阵化为一个正交矩阵Q和一个上三角矩阵R左乘的形式。
从QR算法的构造过程可以看到算法的主要计算量出现在QR分解上,如果直接对矩阵A用QR方法求全部特征值,那麽涉及的计算量是很大的,因此应该先对A作预处理。
应用中常先对做正交相似变换将其化为上Hessenberg矩阵H,然后再对H采用QR方法,可以大大减少计算量,这里Hessenberg矩阵也称为拟三角矩阵,它的非零元素比三角矩阵多了一条次对角线。
采用此程序计算题目所给矩阵的特征值时,迭代次数k=20,可见采用带双步位移的QR算法,可以减少迭代次数加速收敛,避免复数运算,在计算上是比较经济的。
6参考文献
【1】冯康等编:
数值计算方法,国防工业出版社,1978年。
【2】颜庆津编:
数值分析(第三版),北京航空航天大学出版社,2006年。
【3】谭浩强编:
C程序设计(第二版),清华大学出版社,2004年。
7附录
附录A全部源程序(包括主程序和每个子程序的功能注释)
///////////////////带双步位移的QR算法求矩阵特征值和特征向量程序//////////////////////
//声明:
本程序为北航2008级研究生《数值分析A》计算实习题第一题。
#include
#include
voidinput(intm,doublen[][10]) //数据输入子程序(将aij的数据输入二维数组)
{
inti,j;
for(i=1;i {
for(j=1;j {
if(i!
=j)
n[i-1][j-1]=sin(0.5*i+0.2*j);
else
n[i-1][j-1]=1.5*cos(i+1.2*j);
}
}
}
intsign(doublea) //定义符号函数的子程序
{
intb;
if(a>0)
b=1;
else
b=-1;
returnb;
}
voidmul(intn1,intn2,intn3,doublea[][10],doublem1[][10],doublem2[][10])//定义矩阵乘法的子程序
{
inti,j,l;
for(i=0;i {
for(j=0;j {
a[i][j]=0;
for(l=0;l {
a[i][j]=a[i][j]+m1[i][l]*m2[l][j];
}
}
}
}
voidequ(doubleA[][10],doubleB[][10])//定义矩阵赋值的子程序
{
inti,j;
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
A[i][j]=B[i][j];
}
}
}
voidsub(intn1,intn2,doublem,doublea[][10],doublem1[][10],doublem2[][10])//定义矩阵减法的子程序
{
inti,j;
for(i=0;i {
for(j=0;j {
a[i][j]=0;
a[i][j]=m1[i][j]-m*m2[i][j];
}
}
}
voidadd(intn1,intn2,doublem,doublea[][10],doublem1[][10],doublem2[][10])//定义矩阵加法的子程序
{
inti,j;
for(i=0;i {
for(j=0;j {
a[i][j]=0;
a[i][j]=m1[i][j]+m*m2[i][j];
}
}
}
voidfunction(doublex[][2],doubleD[][2]) //求解二次方程X*X+a*X+b=0的子程序
{
do