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

《数值分析A》 第一次数值分析大作业 姓名： 王 硕 学号： SY1415309 一．算法设计方案： 1．由于λ1＜λ2＜…＜λ501，所以在以所有特征值建立的数轴上，λ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的和。 6．A的谱范数条件数等于按模最大的特征值的绝对值与按模最小的特征值的绝对值之比，按模最大的特征值与按模最小的特征值已分别在前面求出。 7．A的行列式的值就是矩阵A进行LU分解后U的对角线元素的乘积。 二．源程序： 该计算程序使用FORTRAN语言编写 module linear1 !生成模块，便于主程序中调用 implicit none contains subroutine lu(matrix,us,b,x) !这个子程序是用来对矩阵进行LU分解的，并且可以求解线性方程AX=B real(kind=8) matrix(:,:),b(:),us(:,:),x(:) integer n real(kind=8),allocatable :: l(:,:),y(:) real(kind=8) a,c,e,f integer i,j,k,t n=size(matrix,1) allocate(l(n,n),y(n)) a=0.0 c=0.0 do j=1,n us(1,j)=matrix(1,j) end do do i=2,n l(i,1)=matrix(i,1)/us(1,1) end do do k=2,n-1 do j=k,n do i=k+1,n do t=1,k-1 a=a+l(k,t)*us(t,j) c=c+l(i,t)*us(t,k) end do us(k,j)=matrix(k,j)-a l(i,k)=(matrix(i,k)-c)/us(k,k) !矩阵A的LU分解结束 a=0.0 c=0.0 end do end do end do a=0.0 do t=1,n-1 a=a+l(n,t)*us(t,n) end do us(n,n)=matrix(n,n)-a y(1)=b(1) e=0.0 f=0.0 do i=2,n do t=1,i-1 e=e+l(i,t)*y(t) end do y(i)=b(i)-e e=0.0 end do x(n)=y(n)/us(n,n) do i=n-1,1,-1 do t=i+1,n f=f+us(i,t)*x(t) end do x(i)=(y(i)-f)/us(i,i) !求得解向量X(i) f=0.0 end do end subroutine lu real function uu(aa,bb) !这个子程序是用来做两个向量啊a,b的乘积的 real(kind=8) aa(:),bb(:) integer n,i n=size(aa) uu=0.0 do i=1,n uu=uu+aa(i)*bb(i) end do end function subroutine au(matrix,u,c) !这个子程序是用来计算矩阵和向量相乘的，最后返回一个向量值 real(kind=8) matrix(:,:),u(:),c(:) integer n,i,k n=size(matrix,1) do i=1,n c(i)=0.0 end do do k=1,n do i=1,n c(k)=c(k)+matrix(k,i)*u(i) end do end do end subroutine au