北航数值分析大作业2Word文档格式.docx

1、这个子程序是用来做两个向量啊a,b的乘积的real（kind=8） aa（:）,bb（:integer n,in=size（aa）uu=0.0uu=uu+aa（i）*bb（i）end functionsubroutine uptr（ma,ma1） !这个子程序是用来对矩阵A进行拟上三角化的real（kind=8） ma（:）,ma1（:integer nreal（kind=8）,allocatable : u（:）,p（:）,q（:）,w（:）,wu（:）,up（:）,at（:real（kind=8） t,d,c0,h,heinteger i,j,rn=size（ma,1）allocate（u

2、（n）,p（n）,q（n）,w（n）,wu（n,n）,up（n,n）,at（n,n）,c（n）do j=1,nma1（i,j）=0.0end do !ma1的零初始化do r=1,n-2he=0.0do i=r+1,nhe=he+ma（i,r）*2d=sqrt（he）if（ma（r+1,r）0） thenc0=0.0-delsec0=dend ifdo i=1,rma1（i,r）=ma（i,r）ma1（r+1,r）=c0do i=r+2,nif（he=0） thencycleh=c0*2-c0*ma（r+1,r）u（i）=0.0u（r+1）=ma（r+1,r）-c0u（i）=ma（i,r）at

3、（i,j）=ma（j,i）call au（at,u,c,n）p（i）=c（i）/hcall au（ma,u,c,n）q（i）=c（i）/ht=uu（p,u）/hw（i）=q（i）-t*u（i）ma（i,j）=ma（i,j）-w（i）*u（j）-u（i）*p（j）ma1（i,n-1）=ma（i,n-1）ma1（i,n）=ma（i,n）end subroutine uptrsubroutine two（d,s1,s2,a1,a2） !这个子程序是用来求二阶子阵的两个特征值的real（kind=8） d（:complex（kind=8） s1,s2real（kind=8） a1,a2,b1,b2n=

4、size（d,1） if（d（n-1,n-1）+d（n,n）*2-4*（d（n-1,n-1）*d（n,n）-d（n-1,n）*d（n,n-1）h=c0*2-c0*b（r,r）u（r）=b（r,r）-c0u（i）=b（i,r）bt（i,j）=b（j,i）c1t（i,j）=c1（j,i）call au（bt,u,c,m）v（i）=c（i）/hb（i,j）=b（i,j）-u（i）*v（j）call au（c1t,u,c,m）call au（c1,u,c,m）c1（i,j）=c1（i,j）-w（i）*u（j）-u（i）*p（j）ak1（i,j）=c1（i,j）end subroutine mkaksu

5、broutine qr（ma,ma1,q1,rq） !这个子程序是用来对任一方阵进行QR分解的，并且返回Q,R的值）,q1（:）,rq（:real（kind=8） d,c0,h,heallocate（u（n）,p（n）,w（n）,at（n,n）,c（n）ma1（i,j）=ma（i,j）q1（i,j）=0.0else q1（i,j）=1.0end if do r=1,n-1he=he+ma1（i,r）*2do i=r,nif（ma1（r,r）h=c0*2-c0*ma1（r,r）u（r）=ma1（r,r）-c0u（r）=ma1（i,r）at（i,j）=ma1（j,i）call au（q1,u,c,

6、n）w（i）=c（i）q1（i,j）=q1（i,j）-w（i）*u（j）/hma1（i,j）=ma1（i,j）-u（i）*p（j）rq（i,j）=0.0do r=1,nrq（i,j）=rq（i,j）+ma1（i,r）*q1（r,j）end subroutine qrsubroutine opposite（matrix,u0,y,bo） !这个子程序是用反幂法求出矩阵的模最小的特征值 real（kind=8） matrix（:）,u0（:）,y（:）,bo integer n,i real（kind=8）,allocatable :）,us（:）,x（: real（kind=8） c,nn bo

7、=100 n=size（matrix,1） allocate（u（n）,us（n,n）,x（n） do i=1,n u（i）=u0（i） end do do while（10） c=bo nn=sqrt（uu（u,u） y（i）=u（i）/nn call lu（matrix,us,y,x） !调用对矩阵进行LU分解的子程序 u（i）=x（i） end do bo=uu（y,u） if （abs（bo-c）/abs（bo）0.0）200 if（abs（a（m,m-1）1）then goto 200d（1,1）=a（m-1,m-1）d（1,2）=a（m-1,m）d（2,1）=a（m,m-1）d（2,2）=a（m,m）call two（d,s1,s2,a11,a2） !求二阶子阵的两个特征值if（m