数值分析
—计算实习作业一
学 姓 学
院： 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)