北航数值分析大作业—程序1.doc

数值分析A第一次数值分析大作业姓名： 王 硕学号： SY1415309一算法设计方案：1由于12501，所以在以所有特征值建立的数轴上，1、501位于数轴的两端，两者之一必为按模最大。利用幂法，可以求出来按模最大的特征值，但该值可能为1也可能为501；2上步求出按模最大的特征值M后，将原矩阵平移-M,再利用幂法求一次平移后矩阵的按模最大的特征值，即是数轴上另一端点值M。3比较M与M+M的大小，大的为矩阵A的最大特征值，小的为A矩阵的最小的特征值；4利用反幂法，求矩阵A的按模最小的特征值。但是反幂法中要用到线性方程组的求解，而原矩阵A又是带状矩阵，采用LU分解。所以在这之前要定义一个LU分解子程序，将A矩阵分解为单位下三角矩阵L和上三角矩阵U的乘积。5先利用循环求出k从1到39变化的uk的值。对于每一个uk值，都将压缩后的原始矩阵A的第三行减去相应的uk，形成新的矩阵A(k)，然后调用LU分解的子程序，利用反幂法求出与uk最接近的特征值，该特征值等于利用反幂法求出的值与uk的和。6A的谱范数条件数等于按模最大的特征值的绝对值与按模最小的特征值的绝对值之比，按模最大的特征值与按模最小的特征值已分别在前面求出。7A的行列式的值就是矩阵A进行LU分解后U的对角线元素的乘积。二源程序：该计算程序使用FORTRAN语言编写module linear1 !生成模块，便于主程序中调用 implicit nonecontains subroutine lu(matrix,us,b,x) !这个子程序是用来对矩阵进行LU分解的，并且可以求解线性方程AX=Breal(kind=8) matrix(:,:),b(:),us(:,:),x(:)integer nreal(kind=8),allocatable : l(:,:),y(:)real(kind=8) a,c,e,finteger i,j,k,tn=size(matrix,1)allocate(l(n,n),y(n)a=0.0c=0.0do j=1,nus(1,j)=matrix(1,j)end do do i=2,nl(i,1)=matrix(i,1)/us(1,1)end do do k=2,n-1do j=k,ndo i=k+1,ndo t=1,k-1a=a+l(k,t)*us(t,j)c=c+l(i,t)*us(t,k)end do us(k,j)=matrix(k,j)-al(i,k)=(matrix(i,k)-c)/us(k,k) !矩阵A的LU分解结束a=0.0c=0.0end doend do end doa=0.0do t=1,n-1a=a+l(n,t)*us(t,n)end do us(n,n)=matrix(n,n)-ay(1)=b(1)e=0.0f=0.0do i=2,ndo t=1,i-1e=e+l(i,t)*y(t)end do y(i)=b(i)-ee=0.0 end dox(n)=y(n)/us(n,n)do i=n-1,1,-1do t=i+1,nf=f+us(i,t)*x(t)end dox(i)=(y(i)-f)/us(i,i) !求得解向量X(i)f=0.0end doend subroutine lu real function uu(aa,bb) !这个子程序是用来做两个向量啊a,b的乘积的real(kind=8) aa(:),bb(:)integer n,in=size(aa)uu=0.0do i=1,nuu=uu+aa(i)*bb(i)end do end function subroutine au(matrix,u,c) !这个子程序是用来计算矩阵和向量相乘的，最后返回一个向量值real(kind=8) matrix(:,:),u(:),c(:)integer n,i,kn=size(matrix,1)do i=1,nc(i)=0.0end do do k=1,ndo i=1,nc(k)=c(k)+matrix(k,i)*u(i) end do end doend subroutine au subroutine positive(matrix,u0,bp) !这个子程序是用幂法计算矩阵的模最大的特征值 real(kind=8) matrix(:,:),u0(:),bp integer n,i real(kind=8),allocatable : u(:),y(:),c(:) real(kind=8) co,nn n=size(matrix,1) allocate(u(n),y(n),c(n) bp=0.0 do i=1,n u(i)=u0(i) end do do while(10) co=bp nn=uu(u,u)*0.5 do i=1,n y(i)=u(i)/nn end do call au(matrix,y,c) !调用矩阵和向量相乘的子程序 do i=1,n u(i)=c(i) end do bp=uu(y,u) write(*,(e20.9) bp if (abs(bp-co)/abs(bp)0) c=bo nn=sqrt(uu(u,u) do i=1,n y(i)=u(i)/nn end do call lu(matrix,us,y,x) !调用对矩阵进行LU分解的子程序 do i=1,n u(i)=x(i) end do bo=uu(y,u) write(*,(e20.9) bo if (abs(bo-c)/abs(bo)r1)thent501=r2t1=r1elset501=r1t1=r2end ifopen(unit=30,file=output/thicknessandf.txt) !将程序运行的结果存储在一个TXT的文档里面write(30,(最小的特征值为t1=,e20.9, 最大的特征值为t501=,e20.9) t1,t501write(*,(最小的特征值为t1=,e20.9, 最大的特征值为t501=,e20.9) t1,t501 call opposite(a0,u0,bo) !调用反幂法子程序，求矩阵A按模最小的特征值ts=bo*(0.0-1.0)write(30,(模最小的特征值为,e20.9) tswrite(*,(模最小的特征值为,e20.9) tsdo k=1,39w=t1+k*(t501-t1)/40.0do i=1,mc0(i,i)=a0(i,i)-wend docall opposite(c0,u0,bo) !对于每一个u(k)调用反幂法，求得原点平移后矩阵模最小特征值t=bo*(0.0-1.0)+w !原点平移后矩阵模最小特征值加上平移量就可以得到所要求特征值write(30,(与U,i2,=,e20.9,最接近的特征值为,e20.9) k,w,twrite(*,(与U,i2,=,e20.9,最接近的特征值为,e20.9) k,w,tend docond2=abs(tm/ts) !求矩阵的二范数write(30,(矩阵A的条件数cond(A)2=,e20.9) cond2 write(*,(矩阵A的条件数cond(A)2=,e20.9) cond2 call lu(a0,us,u0,x) det=1.0do i=1,mdet=det*us(i,i) !求矩阵的特征值end do write(30,(矩阵A的行列式DET(A)=,e20.9) detwrite(*,(矩阵A的行列式DET(A)=,e20.9) detstopend program !主程序结束三输出结果：最小的特征值为t1= -0.107000933E+02 最大的特征值为t501= 0.972455120E+01模最小的特征值为 -0.555854581E-02与U 1= -0.101894772E+02最接近的特征值为 -0.101829363E+02与U 2= -0.967886105E+01最接近的特征值为 -0.958562413E+01与U 3= -0.916824493E+01最接近的特征值为 -0.917267607E+01与U 4= -0.865762882E+01最接近的特征值为 -0.865228622E+01与U 5= -0.814701271E+01最接近的特征值为 -0.809346541E+01与U 6= -0.763639660E+01最接近的特征值为 -0.765940833E+01与U 7= -0.712578049E+01最接近的特征值为 -0.711968532E+01与U 8= -0.661516438E+01最接近的特征值为 -0.661175517E+01与U 9= -0.610454826E+01最接近的特征值为 -0.606611000E+01与U10= -0.559393215E+01最接近的特征值为 -0.558510429E+01与U11= -0.508331604E+01最接近的特征值为 -0.511407297E+01与U12= -0.457269993E+01最接近的特征值为 -0.457887511E+01与U13= -0.406208382E+01最接近的特征值为 -0.409651422E+01与U14= -0.355146770E+01最接近的特征值为 -0.355421889E+01与U15= -0.304085159E+01最接近的特征值为 -0.304108999E+01与U16= -0.253023548E+01最接近的特征值为 -0.253397675E+01与U17= -0.201961937E+01最接近的特征值为 -0.200322594E+01与U18= -0.150900326E+01最接近的特征值为 -0.150355918E+01与U19= -0.998387146E+00最接近的特征值为 -0.993558698E+00与U20= -0.487771034E+00最接近的特征值为 -0.487042661E+00与U21= 0.228450775E-01最接近的特征值为 0.223074975E-01与U22= 0.533461189E+00最接近的特征值为 0.532417509E+00与U23= 0.104407730E+01最接近的特征值为 0.105290380E+01与U24= 0.155469341E+01最接近的特征值为 0.158944712E+01与U25= 0.206530952E+01最接近的特征值为 0.206034137E+01与U26= 0.257592564E+01最接近的特征值为 0.255807431E+01与U27= 0.308654175E+01最接近的特征值为 0.308024567E+01与U28= 0.359715786E+01最接近的特征值为 0.361362115E+01与U29= 0.410777397E+01最接近的特征值为 0.409136867E+01与U30= 0.461839008E+01最接近的特征值为 0.460302209E+01与U31= 0.512900620E+01最接近的特征值为 0.513292891E+01与U32= 0.563962231E+01最接近的特征值为 0.559490416E+01与U33= 0.615023842E+01最接近的特征值为 0.608092462E+01与U34= 0.666085453E+01最接近的特征值为 0.668036226E+01与U35= 0.717147064E+01最接近的特征值为 0.729389173E+01与U36= 0.768208675E+01最接近的特征值为 0.771711848E+01与U37= 0.819270287E+01最接近的特征值为 0.822521984E+01与U38= 0.870331898E+01最接近的特征值为 0.864865720E+01与U39= 0.921393509E+01最接近的特征值为 0.925419879E+01矩阵A的条