北航数值分析大作业一
北航数值分析大作业 第一题 幂法与反幂法
数 值 分 析（B ） 大 作 业（一）姓名： 学号： 电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。
1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。
b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。
c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。
②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。
在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。
③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。
求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。
2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。
北航数值分析实验报告
北航数值分析实验报告篇一：北航数值分析报告第一大题《数值分析》计算实习报告第一大题学号：DY1305姓名：指导老师：一、题目要求已知501*501阶的带状矩阵A，其特征值满足?1?2...?501。
试求：1、?1，?501和?s的值；2、A的与数?k??1?k?501??140最接近的特征值?ik（k=1，2，...，39）；3、A的（谱范数）条件数c nd(A)2和行列式de tA。
二、算法设计方案题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵A1如下所示。
?0?0?A1??a1??b??c0b a2bcc bb c............c bb ccb a500b0a 3...a499c?b??a501??0?0??由矩阵A的特征值满足的条件可知?1与?501之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为?501；否则为?1。
在求得?1与?501其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。
用反幂法求得的特征值必为?s。
由条件数的性质可得，c nd(A)2为模最大的特征值与模最小的特征值之比的模，因此，求出?1，?501和?s的值后，则可以求得c nd(A)2。
北航研究生数值分析编程大作业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λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R 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 λ。
北航数值分析计算实习1
《数值分析》计算实习题目110091013 劳云杰一、算法设计方案根据提示的算法，首先使用幂法求出按模最大的特征值λt1，再根据已求出的λt1用带原点平移的幂法求出另一个特征值λt2，比较两个λ的大小，根据已知条件，可以得出λ1和λ501.至于λs，由于是按模最小的特征值，使用反幂法求之，由于反幂法需要解线性方程组，故对矩阵进行Doolittle分解。
再通过带原点平移的反幂法求跟矩阵的与数最接近的特征值。
对非奇异的矩阵A，根据条件数定义，取λt1/λs的绝对值，两个特征值在之前步骤中均以求得。
由于对矩阵进行了Doolittle分解，所以矩阵的行列式det A可由分解得出的上三角阵U 的对角线上元素相乘求得。
为了使A的所有零元素都不存储，使用书本25页的压缩存储法对A进行存储，在计算时通过函数在数组C中检索A中元素即可。
由于A是501*501矩阵，C应取为5*501矩阵。
由于数据不大，为了方便起见，在程序中取502*502矩阵或者502向量，C也取为6*502矩阵。
程序编写参考《数值分析》颜庆津著和[C数值算法].(美国)W ILLIAM.H．P RESS.扫描版。
二、全部源程序#include <stdio.h>#include <math.h>#define XS 1.0e-12//精度水平void fz_a();//对矩阵A赋值double js(int,int);//在压缩矩阵中检索A的元素double mf(double);//幂法double fmf(double);//反幂法int lu(double);//Doolittle分解int jfc(double[],double[]);//解方程int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角阵double (*l)[502]=new double[502][502];//单位下三角阵double a[6][502];//压缩存储矩阵int max(int x,int y)//比大小函数×2{ return (x>y?x:y);}int min(int x,int y)//精度关系，比较下标用{ return (x<y?x:y);}int main(){printf("请耐心等待,先看看中间过程吧～\n");int i,k;double ldt1,ldt2,ld1,ld501,lds,mu[40],det;double ld[40];fz_a();//对A赋值ldt1=mf(0);//幂法求模最大的特征值ldt2=mf(ldt1);//以第一次求得的特征值进行平移ld1=ldt1<ldt2?ldt1:ldt2;//大的就是λ501ld501=ldt1<ldt2?ldt2:ldt1;lu(0);lds=fmf(0);//反幂法求λsdet=1;//初始化行列式for(i=1;i<=501;i++)det=det*u[i][i];//用U的对角元素求行列式for(k=1;k<=39;k++){mu[k]=ld1+k*(ld501-ld1)/40;//与数lu(mu[k]);ld[k]=fmf(mu[k]);}printf("\n 列出结果\n");printf("λ1=%1.12e λ501=%1.12e\n",ld1,ld501);printf("λs=%1.12e \n",lds);printf("cond(A)=%1.12e \n",fabs(ldt1/lds));printf("detA=%1.12e \n",det);for(k=1;k<=39;k++)//列出跟与数最接近特征值{printf("λi%d=%1.12e\t",k,ld[k]);if(k%2==0)printf("\n");}//界面友好性delete []u;delete []l;getchar();return 0;}void fz_a()//对A赋值{int i;for(i=3;i<=501;i++)a[1][i]=a[5][502-i]=-0.064;//原A矩阵的cfor(i=2;i<=501;i++)a[2][i]=a[4][502-i]=0.16;//原A矩阵的bfor(i=1;i<=501;i++)a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);//原对角线元素}double js(int i,int j)//对压缩矩阵检索A的元素{if(abs(i-j)<=2)return a[i-j+3][j];else return 0;}double mf(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;//用幂法的第一种迭代方法for(i=1;i<=501;i++) //用到了2-范数u[i]=1,y[i]=0;for(int k=1;k<=10000;k++)//对迭代次数进行限制{yita=0;for(i=1;i<=501;i++)yita=sqrt(yita*yita+u[i]*u[i]);for(i=1;i<=501;i++)y[i]=u[i]/yita;for(x1=1;x1<=501;x1++){u[x1]=0;for(int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(js(x1,x2)-offset):js(x1,x2))*y[x2];}prebeta=beta;beta=0;for(i=1;i<=501;i++)beta=beta+y[i]*u[i];if(fabs((prebeta-beta)/beta)<=XS){printf("offset=%f lb=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};}//满足误差条件后，迭代终止，并输出平移量，误差和迭代次数return(beta+offset);//加上平移量，方便比较}double fmf(double offset)//反幂法{ int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for(i=1;i<=501;i++)u[i]=1,y[i]=0; //相关量初始化for(int k=1;k<=10000;k++)//限制迭代次数{yita=0;for(i=1;i<=501;i++)yita=sqrt(yita*yita+u[i]*u[i]);for(i=1;i<=501;i++)y[i]=u[i]/yita;jfc(u,y);prebeta=beta;beta=0;for(i=1;i<=501;i++)beta=beta+y[i]*u[i];beta=1/beta;if(fabs((prebeta-beta)/beta)<=XS){printf("offset=%f lb=%f err=%ek=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};}//满足误差条件后，迭代终止，并输出平移量，误差和迭代次数return(beta+offset);}int lu(double offset)//Doolittle分解{int i,j,k,t;double sum;//中间量for(k=1;k<=501;k++)for(j=1;j<=501;j++){u[k][j]=l[k][j]=0;if(k==j)l[k][j]=1;}//对LU矩阵初始化for(k=1;k<=501;k++)//对式（2.12）的程序实现{for(j=k;j<=min(k+2,501);j++){sum=0;for(t=max(1,max(k-2,j-2));t<=(k-1);t++)sum=sum+l[k][t]*u[t][j];//j=k,k+1,……,nu[k][j]=((k==j)?(js(k,j)-offset):js(k,j))-sum;}if(k==501)continue;for(i=k+1;i<=min(k+2,501);i++)//i=k+1,……，n{sum=0;for(t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(js(i,k)-offset):js(i,k))-sum)/u[k][k];}}return 0;}int jfc(double x[],double b[])//解方程{int i,t;double y[502];double sum;y[1]=b[1];for(i=2;i<=501;i++){sum=0;for(t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for(i=500;i>=1;i--){sum=0;for(t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}三、结果λ1=-1.070011361502e+001λ501=9.724634098777e+000λs=-5.557910794230e-003cond(A)=1.925204273902e+003detA=2.772786141752e+118λi1=-1.018293403315e+001 λi2=-9.585707425068e+000 λi3=-9.172672423928e+000λi4=-8.652284007898e+000 λi5=-8.0934********e+000 λi6=-7.659405407692e+000λi7=-7.119684648691e+000 λi8=-6.611764339397e+000 λi9=-6.0661********e+000λi10=-5.585101052628e+000 λi11=-5.114083529812e+000 λi12=-4.578872176865e+000λi13=-4.096470926260e+000 λi14=-3.554211215751e+000 λi15=-3.0410********e+000 λi16=-2.533970311130e+000 λi17=-2.003230769563e+000 λi18=-1.503557611227e+000 λi19=-9.935586060075e -001 λi20=-4.870426738850e -001 λi21=2.231736249575e -002 λi22=5.324174742069e -001 λi23=1.052898962693e+000 λi24=1.589445881881e+000 λi25=2.060330460274e+000 λi26=2.558075597073e+000 λi27=3.080240509307e+000 λi28=3.613620867692e+000 λi29=4.0913********e+000 λi30=4.603035378279e+000 λi31=5.132924283898e+000 λi32=5.594906348083e+000 λi33=6.080933857027e+000 λi34=6.680354092112e+000 λi35=7.293877448127e+000 λi36=7.717111714236e+000 λi37=8.225220014050e+000 λi38=8.648666065193e+000 λi39=9.254200344575e+000四、讨论迭代初始向量的选取对计算结果的影响1.在反幂法中取迭代向量u[1]=1,u[i]=0，i=2,……，501，最后得出的结果中λs=2.668886923785e -002，cond(A)也随之改变成4.009204556274e+0022.在幂法中取迭代向量u[1]=1,u[i]=2,i=2,……,501,最后得出的结果不变。
数值分析大作业一
数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。
再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。
2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。
3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。
4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。
二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。
北航数值分析大作业一
《数值分析B》大作业一SY1103120 朱舜杰一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] .由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是：A的带内元素a ij=C中的元素c i-j+2,j2.求解λ1，λ501，λs①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。
λmin即为λs；如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。
②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max，如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。
3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。
使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。
4.求解A的（谱范数）条件数cond（A）2和行列式d etA。
①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。
二．源程序#include<stdio.h>#include<iostream.h>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip.h>#include<time.h>#define E 1.0e-12 /*定义全局变量相对误差限*/int max2(int a,int b) /*求两个整型数最大值的子程序*/{if(a>b)return a;elsereturn b;}int min2(int a,int b) /*求两个整型数最小值的子程序*/{if(a>b)return b;elsereturn a;}int max3(int a,int b,int c) /*求三整型数最大值的子程序*/{ int t;if(a>b)t=a;else t=b;if(t<c) t=c;return(t);}void assignment(double array[5][501]) /*将矩阵A转存为数组C[5][501]*/{int i,j,k;//所有元素归零for(i=0;i<=4;){for(j=0;j<=500;){array[i][j]=0;j++;}i++;}//第0,4行赋值for(j=2;j<=500;){k=500-j;array[0][j]=-0.064;array[4][k]=-0.064;j++;}//第1,3行赋值for(j=1;j<=500;){k=500-j;array[1][j]=0.16;array[3][k]=0.16;j++;}//第2行赋值for(j=0;j<=500;){ k=j;j++;array[2][k]=(1.64-0.024*j)*sin((double)(0.2*j))-0.64*exp((double)(0.1/j));}}double mifa(double u[501],double array[5][501],double p) /*带原点平移的幂法*/ {int i,j; /* u[501]为初始迭代向量*/double a,b,c=0; /* array[5][501]为矩阵A的转存矩阵*/double y[501]; /*p为平移量*/for(;;){a=0;b=0;/*选用第一种迭代格式*///求ηk-1for(i=0;i<=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求y k-1for(i=0;i<=500;i++){y[i]=u[i]/a;}//求u kfor(i=0;i<=500;i++){u[i]=0;for(j=max2(i-2,0);j<=min2(i+2,500);j++){u[i]+=array[i-j+2][j]*y[j];}u[i]=u[i]-p*y[i]; /*引入平移量*/}//求βkfor(i=0;i<=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)<=E) /*达到精度水平，迭代终止*/break;c=b;}return (b+p); /*直接返回A的特征值*/}void chuzhi(double a[]) /*用随机数为初始迭代向量赋值*/ {int i;srand((int)time(0));for(i=0;i<=500;i++){a[i]=(10.0*rand()/RAND_MAX); /*生成0~10的随机数*/}}void chuzhi2(double a[],int j) /*令初始迭代向量为e i*/{int i;for(i=0;i<=500;i++){a[i]=0;}a[j]=1;}void LU(double array[5][501]) /*对矩阵A进行Doolittle分解*/{ /*矩阵A转存在C[5][501]中*/int j,k,t; /*分解结果L，U分别存在C[5][501]的上半部与下半部*/ for(k=0;k<=500;k++){for(j=k;j<=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t<=(k-1);t++){array[k-j+2][j]-=array[k-t+2][t]*array[t-j+2][j];}}if(k<500)for(j=k+1;j<=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t<=(k-1);t++){array[j-k+2][k]-=array[j-t+2][t]*array[t-k+2][k];}array[j-k+2][k]=array[j-k+2][k]/array[2][k];}}}double fmifa(double u[501],double array[5][501],double p){ /*带原点平移的反幂法*/ int i,j;double a,b,c=0;double y[501];//引入平移量for(i=0;i<=500;i++){array[2][i]-=p;}//先将矩阵Doolittle分解LU(array);for(;;){a=0;b=0;//求ηk-1for(i=0;i<=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求y k-1for(i=0;i<=500;i++){y[i]=u[i]/a;}//回带过程，求解u kfor(i=0;i<=500;i++){u[i]=y[i];}for(i=1;i<=500;i++){for(j=max2(0,(i-2));j<=(i-1);j++){u[i]-=array[i-j+2][j]*u[j];}}u[500]=u[500]/array[2][500];for(i=499;i>=0;i--){for(j=i+1;j<=min2((i+2),500);j++){u[i]-=array[i-j+2][j]*u[j];}u[i]=u[i]/array[2][i];}//求βkfor(i=0;i<=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)<=E) /*达到精度要求，迭代终止*/break;c=b;}return (p+(1/b)); /*直接返回距离原点P最接近的A的特征值*/ }//主函数main(){ int i;double d1,d501,ds,d,a;double u[501];double MatrixC[5][501];printf(" 《数值分析》计算实习题目第一题\n");printf(" SY1103120 朱舜杰\n");//将矩阵A转存为MatrixCassignment(MatrixC);//用带原点平移的幂法求解λ1，λ501chuzhi(u);d=mifa(u,MatrixC,0);chuzhi(u);a=mifa(u,MatrixC,d);if(d<0){d1=d;d501=a;}else{d501=d;d1=a;}printf("λ1=%.12e\n",d1);printf("λ501=%.12e\n",d501);//用反幂法求λschuzhi(u);ds=fmifa(u,MatrixC,0);printf("λs=%.12e\n",ds);//用带原点平移的反幂法求λikfor(i=1;i<=39;i++){a=d1+(i*(d501-d1))/40;assignment(MatrixC);chuzhi(u);d=fmifa(u,MatrixC,a);printf("与μ%02d=%+.12e最接近的特征值λi%02d=%+.12e\n",i,a,i,d);}//求A的条件数d=fabs((d1/ds));printf("A的（谱范数）条件数cond<A>2=%.12e\n",d);//求detAassignment(MatrixC);LU(MatrixC);a=1;for(i=0;i<=500;i++){a*=MatrixC[2][i];}printf("行列式detA=%.12e\n",a);//测试不同迭代初始向量对λ1计算结果的影响。
北航数值分析大作业一
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有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 ii 。
矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。
A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
A 的(谱范数)条件数2)A (cond 和行列式detA 。
一 方案设计1 求1λ，501λ和s λ的值。
s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。
可使用反幂法求得。
1λ，501λ分别为最大特征值及最小特征值。
可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。
使用位移的方式求得另一特征值即可。
2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。
题目可看成求以k μ为偏移量后，按模最小的特征值。
即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。
3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。
矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。
detA 可由LU 分解得到。
因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。
二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。
北航数值分析计算实习第一题编程
i − t + s +1,t t − k + s +1, k t = max(1,i − r ,k − s )
∑c
c
) / cs +1, k
[i = k + 1, k + 2,⋯ , min( k + r , n); k < n]
（2） 求解 Ly = b,Ux = y (数组 b 先是存放原方程右端向量，后来存放中间向量 y)
0 b a2
b c
c b a3 b c
⋯ ⋯ ⋯ ⋯ ⋯
c b a499 b c
c b a500 b 0
c ⎤ b ⎥ ⎥ a501 ⎥ ⎥ 0 ⎥ 0 ⎥ ⎦
在数组 C 中检索矩阵 A 的带内元素 aij 的方法是： A 的带内元素 aij =C 中的元素 ci − j + s +1, j
2
数值分析计算实习题目一
i −1
bi := bi −
பைடு நூலகம்
i − t + s +1,t t t = max(1,i − r )
∑c
b
(i = 2,3,⋯ , n)
xn := bn / cs +1, n
min( i + s )
xi := (bi −
t = i +1
∑c
i −t + s +1,t t
x ) / cs +1,i
(i = n − 1, n − 2,⋯ ,1)
3、Doolittle 分解求解 n 元带状线性方程组（doolittle（）函数）
按照上述对带状矩阵 A 的存储方法和元素 aij 的检索方法，并且把三角分解的结果 ukj 和 lik 分 别存放在 akj 和 aik 原先的存储单元内，那么用 Doolittle 分解法求解 n 元带状线性方程组的算法 可重新表述如下（其中“:=”表示赋值） ： （1） 作分解 A = LU 。 对于 k=1，2， ……，n 执行
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
—计算实习作业一
学 姓 学
院： xxx 名： xxx 号： xxx
2014-11-16
1. 设计方案
观察矩阵 A，结构为带状，且与主对角线相邻的两个带的值 b 和 c 都是常数。从 而可以用带原点平移的幂法或反幂法计算 λ1 和 λ501。 所以算法的设计方案如下： 1.求按模最大的特征值，并记为 lamda1； 2.平移矩阵得到 A’=A-lamda1I，再次用幂法，这次求出的 A’的按模大的特征值 lamda2 就是与步骤 1 求出的特征值相差最大的特征值。即两者一个为最大的特 征值，另一个为最小的特征值。 3.根据 lamda1 和 lamda2 的正负性，直接确定 λ1，和 λ501。 4.对原矩阵 A 用反幂法，求出其按模最小的特征值，记为 lamdas，此即 λs。
end do gama=sqrt(sum) sum=0.0 do i=1,n ue(i)=u(i)/gama end do f=ue if(g==1)then do k=1,n do j=k,min(k+s,n) do t=max(1,k-r,j-s),k-1 sum=sum+a(k-t+s+1,t)*a(t-j+s+1,j) end do a(k-j+s+1,j)=a(k-j+s+1,j)-sum sum=0.0 end do if(abs(a(s+1,k)-0.0)<=1e-6)then write(*,*)"算法失效" exit end if!算法失效的判断 do i=k+1,min(k+r,n) do t=max(1,i-r,k-s),k-1 sum=sum+a(i-t+s+1,t)*a(t-k+s+1,k) end do a(i-k+s+1,k)=(a(i-k+s+1,k)-sum)/a(s+1,k) sum=0.0 end do end do end if !-------------求解Ly=b,Ux=y方程,求出u-----------do i=2,n do t=max(1,i-r),i-1 sum=sum+a(i-t+s+1,t)*f(t) end do f(i)=f(i)-sum sum=0.0 end do u(n)=f(n)/a(s+1,n) do i=n-1,1,-1 do t=i+1,min(i+s,n) sum=sum+a(i-t+s+1,t)*u(t) end do u(i)=(f(i)-sum)/a(s+1,i)
deta=sum sum=1.0 !-----------------输出所有结果----------------------------------open(11,file='所有结果.txt') write(11,'(A,e20.12)')'最小特征值：',lamda(1) write(11,'(A,e20.12)')'最大特征值：',lamda(n) write(11,'(A,e20.12)')'按模最小特征值：',lamdas write(11,'(A,e20.12)')'按模最大特征值：',lamdam1 write(11,'(A,e20.12)')"与数mu最接近的特征值：" do i=1,39 write(11,'(I2,e20.12)')i,lamdai(i) end do write(11,'(A,e20.12)')'条件数cond（A）：',conda write(11,'(A,e20.12)')'行列式的值detA：',deta close(11) pause end program !*************************************************** !---------------幂法求按模最大特征值-------------subroutine mimethod(a,lamda) use array implicit none real(8),intent(in)::a(m,n) real(8),intent(out)::lamda!a为传递进来的矩阵，lamda为传递进来的特征值 integer::k,i,j,t real(8)::sum=0.0D0 real(8)::u(n),ue(n)!u单位化变为ue real(8)::belta=0.0D0,gama,belta0!belta为特征值的迭代近似值，gama为u的 模,belta0在做误差选择中用到 real(8),parameter::epsu=1e-12 k=1 u=1.0D0 u(1)=1.0D0 do while (.true.) write(*,*)'幂法迭代次数K=',k k=k+1 do i=1,n sum=sum+u(i)*u(i) end do gama=sqrt(sum) sum=0.0D0
cond ( A2 )
6.根据公式 其中 λmax 和 λmin 分别是矩阵 A 的按模最大和按 模最小的特征值。将 lamda1 和 lamdas 分别代入 λmax 和 λmin，就可以求出 cond(A) 。 7.对矩阵 A 做三角分解，A=LU。则有
det A u (i, i )
任取非零向量u 0 R n T k 1 u k 1u k 1 y k 1 u k 1 / k 1 Au k y k 1 T k y k 1u k (k 1,2, )
在反幂法的求解过程中，每迭代一次都要求满足解线性方程组 Auk=yk-1。本题 中矩阵 A 上半带宽为 2， 下半带宽也为 2 。 故选择采用三角分解法求解方程组： 先将原矩阵改写成 5 行 501 列的矩阵 C（不存储 A 的 0 元素） A 的带内元素 aij=c 中的元素 ci-j+3。 再对 C 矩阵做 LU 分解。对于 k=1,2，…，n，执行
c(s+2,i-1)=b end do do i=3,n c(s+3,i-2)=cc end do do i=1,n-1 c(s,i+1)=b end do do i=1,n-2 c(s-1,i+2)=cc end do!输入系数矩阵C do i=1,m write(10,*)(c(i,j),j=1,n) end do close(10) !--------------第一问：求最小最大特征值和按模最小特征值----------call mimethod(c,lamdam1)!lamda1是矩阵A的按模最大的特征值 tempc=c do i=1,n tempc(s+1,i)=c(s+1,i)-lamdam1 end do call mimethod(tempc,lamdam2) lamda(1)=min(lamdam1,lamdam1+lamdam2) lamda(n)=max(lamdam1,lamdam1+lamdam2) tempc=c call antimimethod(tempc,lamdas) tempc=c!反幂法会使tempc发生变化 !-----------------第二问：求中间39个特征值------------------------do k=1,39 mu(k)=lamda(1)+k*(lamda(n)-lamda(1))/40.0D0 end do do k=1,39 do i=1,n tempc(s+1,i)=c(s+1,i)-mu(k) end do call antimimethod(tempc,lamdai(k)) lamdai(k)=lamdai(k)+mu(k) tempc=c end do !-----------------第三问：求A的条件数和行列式---------------------conda=abs(lamdam1/lamdas) do i=1,n sum=sum*ukk(i) end do
do i=1,n ue(i)=u(i)/gama end do do i=1,n do t=max(1,i-r),min(i+s,n) sum=sum+a(i-t+s+1,t)*ue(t) end do u(i)=sum sum=0.0 end do do i=1,n sum=sum+ue(i)*u(i) end do belta0=belta belta=sum sum=0.0 if (abs(belta-belta0)/abs(belta)<=epsu)then exit end if end do lamda=belta return end subroutine mimethod !**********************反幂法求按模最小特征值**************** subroutine antimimethod(a,lamda) use array implicit none real(8)::a(m,n) real(8),intent(out)::lamda!a为传递进来的矩阵，lamda为传递进来的特征值 integer::k,i,j,g,t real(8)::sum=0.0 real(8)::u(n),ue(n),f(n)!u单位化变为ue real(8)::belta,gama,belta0!belta为特征值的迭代近似值，gama为u的模,belta0 在做误差选择中用到 real(8),parameter::epsu=1e-12 g=1 belta=0.0D0 u=1.0 do while (.true.) write(*,*)'反幂法迭代次数G=',g do i=1,n sum=sum+u(i)*u(i)