北航数值分析全部三次大作业.pdf
《北航数值分析全部三次大作业.pdf》由会员分享,可在线阅读,更多相关《北航数值分析全部三次大作业.pdf(44页珍藏版)》请在淘文阁 - 分享文档赚钱的网站上搜索。
1、北航数值分析计算实习题目一北航数值分析计算实习题目一一、算法设计:1.要计算矩阵的最大最小特征值,可先通过幂法求得模最大的特征值a,再根据a进行原点平移求出下一特征值b,比较两值的大小,循环下去,最终1和最大特征值501。可以确定最小特征值2.要计算矩阵按模最小的特征值s,可以用反幂法求得,此处需要求解中间过渡向量,采用 LU 分解法解带状方程组,故要设计 Doolite 分解。3.要求解与给定数值接近的特征值,可以对该数做漂移量,新数组特征值倒数的绝对值满足反幂法的要求,故通过反幂法即可求得。4.要计算 cond(A),可由公式 cond(A)=其中max和minmax/min求得,max分
2、别矩阵 A 的模为最大和模为最小的特征值,min可已通过反幂法求得。可以通过幂法求得,而5.要计算 detA,可对 A 作 LU 分解,即 A=LU,则A=LU,其中L 为单位下三角阵,U 为上三角阵,则 detA 即为 U 矩阵主对角元素之积。6.由于矩阵 A 中的零元素较多,为了节省存储空间,可以将矩阵 A 中的元素存储在 5*501 的二维数组中。编译环境:vc+6.0编译函数:幂法,反幂法,Doolite 分解二、全部源程序:#include#include#include#includeint Max(int v1,int v2);int Min(int v1,int v2);voi
3、d zhuancun(double A5501);double mifa(double A5501);void doolite(double A5501,double x501,double b501);double fanmifa(double A5501);double Det(double A5501);void main()double b=0.16,c=-0.064,p,q;int i,j;double A5501;zhuancun(A,b,c);/进行 A 的赋值cout.precision(12);/定义输出精度double lamda1,lamda501,lamdas;doub
4、le k=mifa(A);if(k0)/判断求得最大以及最小的特征值.如果 K0,则它为最大特征值值,/并以它为偏移量再用一次幂法求得新矩阵最大特征值,即为最大 /与最小的特征值的差lamda501=k;for(i=0;i=500;i+)A2i=A2i-k;lamda1=mifa(A)+lamda501;for(i=0;i=500;i+)A2i=A2i+k;else /如果 K=0,则它为最小特征值值,并以它为偏移量再用一次幂法 /求得新矩阵最大特征值,即为最大与最小的特征值的差 lamda1=k;for(i=0;i=500;i+)A2i=A2i-k;lamda501=mifa(A)+lamd
5、a1;for(i=0;i=500;i+)A2i=A2i+k;lamdas=fanmifa(A);FILE*fp=fopen(result.txt,w);fprintf(fp,1=%.12en,lamda1);fprintf(fp,501=%.12en,lamda501);fprintf(fp,s=%.12enn,lamdas);fprintf(fp,t 要求接近的值ttt 实际求得的特征值n);for(i=1;i=39;i+)/反幂法求得与给定值接近的特征值p=lamda1+(i+1)*(lamda501-lamda1)/40;for(j=0;j=500;j+)A2j=A2j-p;q=fanm
6、ifa(A)+p;for(j=0;jv2)?v1:v2);int Min(int v1,int v2)return(v1v2)?v1:v2);/*将矩阵值转存在一个数组里*/void zhuancun(double A5501,double b,double c)int i=0,j=0;Aij=0,Aij+1=0;for(j=2;j=500;j+)Aij=c;i+;j=0;Aij=0;for(j=1;j=500;j+)Aij=b;i+;for(j=0;j=500;j+)Aij=(1.64-0.024*(j+1)*sin(0.2*(j+1)-0.64*exp(0.1/(j+1);i+;for(j
7、=0;j=499;j+)Aij=b;Aij=0;i+;for(j=0;j=498;j+)Aij=c;Aij=0,Aij+1=0;/*幂法求解模最大的特征值*/double mifa(double A5501)int s=2,r=2,m=0,i,j;double b2,b1=0,sum,u501,y501;for(i=0;i=500;i+)ui=1.0;dosum=0;if(m!=0)b1=b2;m+;for(i=0;i=500;i+)sum+=ui*ui;for(i=0;i=500;i+)yi=ui/sqrt(sum);for(i=0;i=500;i+)ui=0;for(j=Max(i-r,0
8、);j=Min(i+s,500);j+)ui=ui+Ai-j+sj*yj;b2=0;for(i=0;i=exp(-12);/设置计算精度return b2;/*带状 DOOLITE 分解,并且求解出方程组的解*/void doolite(double A5501,double x501,double b501)int i,j,k,t,s=2,r=2;double B5501,c501;for(i=0;i=4;i+)for(j=0;j=500;j+)Bij=Aij;for(i=0;i=500;i+)ci=bi;for(k=0;k=500;k+)for(j=k;j=Min(k+s,500);j+)
9、for(t=Max(0,Max(k-r,j-s);t=k-1;t+)Bk-j+sj=Bk-j+sj-Bk-t+st*Bt-j+sj;for(i=k+1;i=Min(k+r,500);i+)for(t=Max(0,Max(i-r,k-s);t=k-1;t+)Bi-k+sk=Bi-k+sk-Bi-t+st*Bt-k+sk;Bi-k+sk=Bi-k+sk/Bsk;for(i=1;i=500;i+)for(t=Max(0,i-r);t=0;i-)xi=ci;for(t=i+1;t=Min(i+s,500);t+)xi=xi-Bi-t+st*xt;xi=xi/Bsi;/*反幂法求解模最大的特征值*/do
10、uble fanmifa(double A5501)int s=2,r=2,m=0,i;double b2,b1=0,sum=0,u501,y501;for(i=0;i=500;i+)ui=1.0;doif(m!=0)b1=b2;m+;sum=0;for(i=0;i=500;i+)sum+=ui*ui;for(i=0;i=500;i+)yi=ui/sqrt(sum);doolite(A,u,y);b2=0;for(i=0;i=fabs(b1)*exp(-12);return 1/b2;/*矩阵 A 的 LU 分解,其中 U 的主线乘积即为矩阵的DET*/double Det(double A5
11、501)int i,j,k,t,s=2,r=2;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+)Ak-j+sj=Ak-j+sj-Ak-t+st*At-j+sj;for(i=k+1;i=Min(k+r,500);i+)for(t=Max(0,Max(i-r,k-s);t=k-1;t+)Ai-k+sk=Ai-k+sk-Ai-t+st*At-k+sk;Ai-k+sk=Ai-k+sk/Ask;double det=1;for(i=0;i=500;i+)det*=Asi;return det;三、运
12、行结果:1=-1.069936345952e+001 501=9.722283648681e+000 s=-5.557989086521e-003要求接近的值实际求得的特征值 1:-9.678281104107e+000 i1:-9.585702058251e+000 2:-9.167739926402e+000 i2:-9.172672423948e+000 3:-8.657198748697e+000 i3:-8.652284007885e+000 4:-8.146657570993e+000 i4:-8.093482780052e+000 5:-7.636116393288e+000 i
13、5:-7.659405420574e+000 6:-7.125575215583e+000 i6:7:-6.615034037878e+000 i7:8:-6.1 04492860173e+000 i8:9:-5.593951682468e+000 i9:10:-5.083410504763e+000 i10:11:-4.572869327058e+000 i11:12:-4.062328149353e+000 i12:13:-3.551786971648e+000 i13:14:-3.041245793943e+000 i14:15:-2.530704616238e+000 i15:16:-
14、2.020163438533e+000 i16:17:-1.509622260828e+000 i17:18:-9.990810831232e-001 i18:19:-4.885399054182e-001 i19:20:2.200127228676e-002 i20:2.231736249587e 21:5.325424499917e-001 22:1.043083627697e+000 i22:1.052898964020e+000 23:1.553624805402e+000 i23:1.589445977158e+000 24:2.064165983107e+000 i24:2.060
15、330427561e+000 25:2.574707160812e+000 i25:2.558075576223e+000 26:3.085248338516e+000 i26:3.080240508465e+000 27:3.595789516221e+000 i27:3.613620874136e+000 28:4.106330693926e+000 i28:4.091378506834e+000 29:4.616871871631e+000 i29:4.603035354280e+000 30:5.127413049336e+000 i30:5 31:5.637954227041e+00
16、0 i31:5.594906275501e+000 32:6.148495404746e+000 i32:6.080933498348e+000 33:6.659036582451e+000 i33:6.680354121496e+000 34:7.169577760156e+000 i34:7.293878467852e+000 35:7.68011 8937861e+000 i35:7.717111851857e+000 36:8.190660115566e+000 i36:8.225220016407e+000 37:8.701201293271e+000 i37:8.648665837
17、870e+000 38:9.211742470976e+000 i38:9.254200347303e+000 39:9.722283648681e+000 i39:9.7246340-7.119684646576e+000-6.611764337314e+000-6.066103126985e+000-5.585101045269e+000-5.114083539196e+000-4.578872177367e+000-4.096473385708e+000-3.554211216942e+000-3.041090018133e+000-2.533970334136e+000-2.00323
18、0401311e+000-1.503557606947e+000-9.935585987809e-001-4.870426734583e-001-002i21:5.324174742068e-001.132924284378e+00099672e+000cond(A)=1.925042185755e+003detA=2.772786141752e+118四、实验分析:实验发现,当初始向量 u0中 0 元素的个数较多时,计算结果与准确值相差较大;而当增加初始向量 u0中 1 元素的个数,发现其结果更加靠近准确值。所以,初始向量的选取会对计算结果产生一定的影响。另外,初始向量选择不当时将导致迭代中
19、 x1 的系数等于零。但是,由于舍入误差的影响,经若干步迭代后,按照基向量展开时,x1 的系数可能不等于零,把这一向量看作初始向量,用幂法继续求向量序列,仍然会得出预期的结果,不过收敛速度较慢。如果收敛很慢,可改换初始向量,所以初始向量的选取还会影响到迭代的次数。北航数值分析计算实习题目二北航数值分析计算实习题目二一算法设计方案一算法设计方案:1.矩阵的 QR 分解。把矩阵A 分解为一个正交矩阵 Q 与一个上三角矩阵 R 的乘积,称为矩阵A 的正三角分解,简称 QR 分解。QR 分解的算法如下:记A1 A,并记Ar aijnn,令Q1 I(n 阶单位矩阵)对于 r=1,2,n-1 执行r(1)
20、若air(i r 1,r 2,.,n)全为零,则令Qr1 Qr,Ar1 Ar转(5);否则转(2)r(2)计算dr(airn(r)ir)2(r)(r)cr sgn(arr)dr(若arr 0,则取cr dr)(r)hr cr2crarr(3)令ur(0,.,0,arrcr,ar1,r,.,anr)R(4)计算(r)(r)(r)Tnr QrurTQr1 Qrrur/hrpr A ur/hrTAr1 ArurprTr(5)继续当此算法执行完后就得到正交矩阵Q Qn和上三角矩阵R An且有A QR。2.矩阵的 拟上三角化。对实矩阵A 的拟上三角化具体算法如下:记A(1)r A,并记A(r)的第 r
21、列到第 n 列的元素为aij(i 1,2,.,n;j r,r 1,.,n)。对于r 1,2,.,n 2执行(r)(1)若air(i r 2,r 3,.,n)全为零,则令A(r1)A(r),转(5);否则转(2)。(2)计算drir1(an(r)ir)2)(r)cr sgn(ar(r1,r)dr(若ar1,r 0,则取cr dr)hr cr2crar(r1,r(3)令ur(0,.,0,ar1,rcr,ar2,r,.,anr)R(4)计算(r)(r)(r)Tnpr A(r)Tur/hrqr A(r)ur/hrTtr prur/hrr qrtrurTTA(r1)A(r)rururpr(5)继续算法执
22、行完后,就得到与原矩阵A 相似的拟上三角矩阵A3.带双步位移的 QR 方法求实矩阵A Rnn(n1)。全部特征值的具体算法如下:nn(1)使用矩阵的拟上三角化的算法把矩阵A R精度水平 0和迭化最大次数 L。(2)记A1 A(n1)化为拟上三角矩阵A(n1);给定(1)aijnn,令k 1,m n。(k)(k)(3)如果am,m1,则得到 A 的一个特征值amm,置m:m 1(降阶),转(4);否则转(5)。(k)(4)如果m 1,则得到 A 的一个特征值a11,转(11);如果m 0,则直接转(11);如果m 1转(3)。(5)求二阶子阵(k)(k)am1,m1am1,mDk(k)(k)am
23、,m1am,m的两个特征值s1和s2,即计算二次方程(k)(k)2(am1,m1 am,m)detDk 0的两个根s1和s2s1和s2,转(11);否则转(7)。(6)如果m 2,则得到 A 的两个特征值(k)(7)如果am1,m2,则得到 A 的两个特征值s1和s2,置m:m 2(降阶),转(4);否则转(8)。(8)如果 k=L,则计算终止,未得到A 的全部特征值;否则转9。(9)记Akaijmm(i i,j m),计算(k)(k)s am1,m1 am,m(k)(k)(k)(k)t am1,m1am,m am,m1am1,m(k)Mk Ak2 sAktI(I是m阶单位矩阵)Mk QkRk
24、(对Mk作QR分解)TAk1 QkAkQk(10)置k:k 1,转(3)。(11)A 的全部特征值已计算完毕,停止计算。4.算法的主过程:(1)根据公式生成原始矩阵A aij1010(2)对矩阵 A 进行拟上三角化(3)用带双步位移的 QR 方法求矩阵的特征值(4)对每个特征值,解出线性方程组(AI)X 0的所有非零解即得 A 的属于 的全部特征向量。二全部源程序:二全部源程序:#include#include#include#define L 2500#define n 11#define e 0.000000000001int i,j,s,p,k,ik;double Ann,qnn,rnn
25、,rqnn,Inn;double Pn,Wn,un,Qn;double dr,cr,hr,ar,tr;int nR,nC;double s1,t,x,tzRn,tzC2n,sum,Mnn,vn;void juzhenA()/生成矩阵 Afor(i=1;i11;i+)for(j=1;j11;j+)if(j!=i)Aij=sin(0.5*i+0.2*j);elseAij=1.5*cos(i+1.2*j);void hess()for(s=1;sn-2;s+)for(ar=0.0,i=s+2;i0)cr=-dr;else cr=dr;hr=cr*cr-cr*As+1s;for(i=1;i=s;i+)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 全部 三次 作业
限制150内