北航数值分析大作业FORTRAN语言版
北航数值分析大作业二(纯原创,高分版)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令 以及最大迭代步数L；
2.若m≤0，则结束计算，已求出A的全部特征值，判断 或 或m≤2是否成立，成立则转3，否则转4；
3.若 ，则得一个特征值 ，m=m-1，降阶；若 ，则计算矩阵：
的特征值得矩阵A的两个特征值，m=m-2，降阶，转2.；
4.若k≤L，成立则令
k=k+1，转2，否则结束计算，为计算出矩阵A的全部特征值；
fortran 语言编程
fortran 语言编程Fortran 语言编程Fortran（Formula Translation）是一种面向科学计算和工程计算的高级编程语言。
它于1957年诞生于IBM，是最早被广泛采用的科学计算语言之一，目前已经发展到第四个版本（Fortran 2018）。
Fortran是一种编译型语言，它通过编写源代码并使用编译器将其转换成机器语言来执行。
本文将详细介绍Fortran语言的基础知识、语法规则和常用的编程技巧，以帮助读者了解和掌握这门强大的科学计算语言。
第一步：安装Fortran编译器要开始编写和运行Fortran程序，首先需要安装Fortran编译器。
有多种Fortran编译器可供选择，其中最常用的是GNU Fortran（gfortran）和Intel Fortran Compiler（ifort）。
可以从官方网站或其他可信的来源获得这些编译器的安装程序，并按照提示进行安装。
第二步：编写并编译Fortran程序在开始编写Fortran程序之前，需要了解Fortran的基本语法规则。
Fortran使用固定格式或自由格式，固定格式的源代码按照列格式排列，每行的前6列被保留用于行号和注释，从第7列开始是可执行代码。
自由格式没有列格式的限制，更加灵活，但在编译阶段需要指定自由格式。
下面是一个简单的Fortran程序示例，用于计算并输出两个数的和：fortranprogram additionimplicit noneinteger :: a, b, sumprint *, "Enter two numbers:"read *, a, bsum = a + bprint *, "The sum is:", sumend program addition将以上代码保存为一个以.f90为后缀名的文件（例如addition.f90），然后使用编译器将其编译成可执行程序。
北航数值分析报告大作业一
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号 ZY1403140学生许阳教师玉泉日期 2014 年 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.求1λ，501λ和s λ的值。
2.求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
3.求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 值。
北航数值分析大作业三
一、题目：关于x, y, t, u, v, w 的下列方程组0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩1、试用数值方法求出f(x, y)在区域 {(,)|00.8,0.5 1.5}D x y x y =≤≤≤≤上的一个近似表达式,0(,)kr s rsr s p x y cx y ==∑要求(,)p x y 一最小的k 值达到以下的精度10202700((,)(,))10i j i j i j f x y p x y σ-===-≤∑∑其中，0.08,0.50.05i j x i y j ==+。
2、计算****(,),(,)i j i j f x y p x y （i = 1, 2, …，8；j = 1, 2,…，5）的值，以观察(,)p x y 逼近(,)f x y 的效果，其中，*i x =0.1i , *j y =0.5+0.2j 。
说明：1、用迭代方法求解非线性方程组时，要求近似解向量()k x 满足()(1)()12||||/||||10k k k x x x --∞∞-≤2、作二元插值时，要使用分片二次代数插值。
3、要由程序自动确定最小的k 值。
4、打印以下内容：●算法的设计方案。
●全部源程序（要求注明主程序和每个子程序的功能）。
●数表：,,i j x y (,)i j f x y （i = 0,1,2,…,10；j = 0,1,2,…,20）。
●选择过程的,k σ值。
●达到精度要求时的,k σ值以及(,)p x y 中的系数rs c （r = 0,1,…,k;s = 0,1,…,k ）。
●数表：**,,i j x y ****(,),(,)i j i j f x y p x y （i = 1, 2, ...,8；j = 1, 2, (5)。
北航数值分析报告大作业第三题(fortran)
“数值分析“计算实习大作业第三题——SY1415215孔维鹏一、计算说明1、将x i=0.08i，y j=0.5+0.05j分别代入方程组（A.3）得到关于t，u，v，w的的方程组，调用离散牛顿迭代子函数求出与x i，y j对应的t i，u j。
2、调用分片二次代数插值子函数在点(t i,u j)处插值得到z(x i,y j)=f(x i,y j)，得到数表(x i,y j,f(x i,y j))。
3、对于k=1,2,3,4⋯，分别调用最小二乘拟合子函数计算系数矩阵c rs及误差σ，直到满足精度，即求得最小的k值及系数矩阵c rs。
4、将x i∗=0.1i，y j∗=0.5+0.2j分别代入方程组（A.3）得到关于t∗，u∗，v∗，w∗的的方程组，调用离散牛顿迭代子函数求出与x i∗，y j∗对应的t i∗，u j∗，调用分片二次代数插值子函数在点(t i∗,u j∗)处插值得到z∗(x i∗,y j∗)=f(x i∗,y j∗)；调用步骤3中求得的系数矩阵c rs求得p(x i∗,y j∗)，打印数表(x i∗,y j∗,f(x i∗,y j∗),p(x i∗,y j∗))。
二、源程序（FORTRAN）PROGRAM SY1415215DIMENSION X(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21),C(6,6) DIMENSION X1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组，得到z=f（t，u）的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDOENDDO!*************分片二次插值得到z=f（x，y）***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDOENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P（x，y）**************N=11M=21WRITE (1,'(" ","K和σ分别为：")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E)WRITE (1,'(I3,2X,E20.13)') K-1,EIF(E<10E-7) EXITENDDOWRITE(1,'(" ")')WRITE(1,'("系数矩阵Crs（按行）为：")')DO I=1,KDO J=1,KWRITE (1,'(E20.13,2X,\)') C(I,J)ENDDOWRITE (1,"('')")WRITE (*,"('')")ENDDODO I=1,8X1(I)=0.1*IENDDODO J=1,5Y1(J)=0.5+0.2*JENDDODO I=1,8DO J=1,5CALL DISNEWTON_NONLINEAR(X1(I),Y1(J),UX1(I,J),TY1(I,J))ENDDOENDDODO I=1,8DO J=1,5CALL INTERPOLATION(UX1(I,J),TY1(I,J),FXY1(I,J))ENDDOENDDOPXY1=0DO I=1,8DO J=1,5DO II=1,KDO JJ=1,KPXY1(I,J)=PXY1(I,J)+C(II,JJ)*(X1(I)**(II-1))*(Y1(J)**(JJ-1)) ENDDOENDDOENDDOENDDOWRITE(1,'(" ")')WRITE(1,'("数表（x,y,f(x,y),p(x,y)):")')WRITE(1,"(2X,'X',6X,'Y',12X,'F(X,Y)',14X,'P(X,Y)')")DO I=1,8DO J=1,5WRITE(1,'(F5.3,2X,F5.3,2X,E20.13,2X,E20.13)') X1(I),Y1(J),FXY1(I,J),PXY1(I,J) ENDDOWRITE (1,"('')")ENDDOCLOSE (1)END!***********用离散牛顿法求解非线性方程组****************SUBROUTINE DISNEWTON_NONLINEAR(X1,Y1,U,T)PARAMETER (N=4)REAL EPS !EPS为迭代精度，M为最大迭代次数DIMENSION X(N),H(N),Y(N),JA(N,N),E(N),XK(N)REAL(8) JA,X,H,Y,E,XK,U,T,V,W,X1,Y1,E1,E2F1(T,U,V,W)=0.5*COS(T)+U+V+W-X1-2.67F2(T,U,V,W)=T+0.5*SIN(U)+V+W-Y1-1.07F3(T,U,V,W)=0.5*T+U+COS(V)+W-X1-3.74F4(T,U,V,W)=T+0.5*U+V+SIN(W)-Y1-0.79EPS=10E-12M=100X=1.0DO K=1,MH=1!计算Y=F（x）Y(1)=F1(X(1),X(2),X(3),X(4))Y(2)=F2(X(1),X(2),X(3),X(4))Y(3)=F3(X(1),X(2),X(3),X(4))Y(4)=F4(X(1),X(2),X(3),X(4))!计算JA（N,N）E=0.0DO I=1,NDO J=1,NDO JJ=1,NIF(JJ==J) THENE(JJ)=X(JJ)+H(JJ)ELSEE(JJ)=X(JJ)ENDIFENDDOIF(I==1) THENJA(I,J)=(F1(E(1),E(2),E(3),E(4))-F1(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==2) THENJA(I,J)=(F2(E(1),E(2),E(3),E(4))-F2(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==3) THENJA(I,J)=(F3(E(1),E(2),E(3),E(4))-F3(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==4) THENJA(I,J)=(F4(E(1),E(2),E(3),E(4))-F4(X(1),X(2),X(3),X(4)))/H(J) ENDIFENDDOENDDO!求解线性方程组CALL GAUSS(JA,XK,-Y,N)!判断精度CALL NORM(XK,N,E1)CALL NORM(X,N,E2)IF(E1/E2<=EPS) THENT=X(1)U=X(2)EXITELSEDO I=1,NX(I)=X(I)+XK(I)ENDDOENDIFENDDORETURNEND!**********列主元高斯消去法求解线性方程组********* SUBROUTINE GAUSS(A,X,B,N)DIMENSION A(N,N),B(N),X(N),T(N,N),TB(N)REAL M(N,N)REAL(8) A,B,X,T!消元过程DO K=1,N-1TA=A(K,K)TL=KDO L=K+1,NIF ((A(L,K)>TA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIFENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J)ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A)DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I))ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标，W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R*************************IF(U<=X(2)+H/2) THENK=2ELSEIF(U>X(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(U<=X(I)+H/2)) THENK=IENDIFENDDOENDIFIF(V<=Y(2)+T/2) THENR=2ELSEIF(V>Y(M-1)-T/2) THENR=M-1ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(V<=Y(J)+T/2)) THENR=JENDIFENDDOENDIFI=KJ=RLK(1)=(U-X(I))*(U-X(I+1))/(X(I-1)-X(I))/(X(I-1)-X(I+1))LK(2)=(U-X(I-1))*(U-X(I+1))/(X(I)-X(I-1))/(X(I)-X(I+1)) LK(3)=(U-X(I))*(U-X(I-1))/(X(I+1)-X(I))/(X(I+1)-X(I-1)) LR(1)=(V-Y(J))*(V-Y(J+1))/(Y(J-1)-Y(J))/(Y(J-1)-Y(J+1)) LR(2)=(V-Y(J-1))*(V-Y(J+1))/(Y(J)-Y(J-1))/(Y(J)-Y(J+1)) LR(3)=(V-Y(J))*(V-Y(J-1))/(Y(J+1)-Y(J))/(Y(J+1)-Y(J-1))W=0DO K=1,3DO R=1,3W=W+LK(K)*LR(R)*Z(J+R-2,I+K-2)ENDDOENDDORETURNEND!*******************最小二乘拟合子函数************** SUBROUTINE LSFITTING(X,Y,Z,A,N,M,P,Q,DT1)INTEGER P,QDIMENSION X(N),Y(M),Z(N,M),A(P,Q)DIMENSION APX(20),APY(20),BX(20),BY(20),U(20,20),V(20,M) DIMENSION T(20),T1(20),T2(20)REAL(8) X,Y,Z,A,DT1DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDOIF(P>N) P=NIF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDOAPX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*GENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,MV(K,J)=V(K,J)+Z(I,J)*GENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDOU(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*GAPY(2)=APY(2)+(Y(I))*G*G ENDDOAPY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*GDO J=1,PU(J,K)=U(J,K)+V(J,I)*GENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDOENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDODO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2)IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K)ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L)ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)):X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.079 0.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.300 1.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.00 1.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769 -0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454 -0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.262 0.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.150 1.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718 -0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.225 0.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.301 0.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+000.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477 -0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.363 0.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.250 1.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739 -0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.199 0.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517 -0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662 -0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674 -0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.223 0.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.494 0.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.450 1.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.168 0.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+000.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.250 1.644 0.511 0.51E+00 0.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.568 0.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.151 0.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.443 0.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为：0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs（按行）为：0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-010.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02数表（x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00 -0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+000.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+01 0.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.700 1.500 -0.53E+00 -0.80E+000.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北航数值分析大作业第二题(fortran)
!计算A（r+1） DO I=1,N DO J=1,N A(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J) ENDDO ENDDO ENDIF ENDDO RETURN END
!***************符号函数子程序*****************！ FUNCTION SGN(X) REAL(8) X IF(X>0) THEN SGN=1 ELSE IF(X<0) THEN SGN=-1 ELSE IF(X==0) THEN SGN=0 ENDIF END
DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数 组,1为实部，为虚部 REAL(8) A,A1,A2,C,Q,R,CM E=1E-12 L=1000 !精度水平 !迭代最大次数
OPEN(1,FILE='数值分析大作业第二题计算结果.TXT') DO I=1,N DO J=1,N IF(I==J) THEN A(I,J)=1.52*COS(I+1.2*J) ELSE A(I,J)=SIN(0.5*I+0.2*J) ENDIF ENDDO ENDDO A1=A WRITE(*,"('矩阵A为：')") WRITE(1,"('矩阵A为：')") DO I=1,N DO J=1,N WRITE(*,"(2X,E20.13,2X,\)") A(I,J) WRITE(1,"(2X,E20.13,2X,\)") A(I,J) ENDDO WRITE(*,"(' ')") WRITE(1,"(' ')") ENDDO !使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1） CALL HESSENBERG(A,N) WRITE(*,"('拟上三角化后矩阵A（n-1）为：')") WRITE(1,"('拟上三角化后矩阵A（n-1）为：')") DO I=1,N DO J=1,N WRITE(*,"(2X,E20.13,2X,\)") A(I,J) WRITE(1,"(2X,E20.13,2X,\)") A(I,J) ENDDO WRITE(*,"('')") WRITE(1,"('')") ENDDO !计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵 A2=A CALL QRD(A2,N,Q,R)
北航数值分析大作业3(学硕)
《数值分析》作业三院系：机械学院学号：SY1307145姓名：龙安林2013年11 月24 日1. 算法设计1) 开始；2) 计算数组[][]0.08,0.050.5,0,1,2,,10;0,1,2,,20x i i y j j i j ==+=⋯=⋯（）； 3) 将点[][],0,1,2,,10;0,1,2,,20x i y j i j =⋯=⋯（）,（）带入非线性方程组： 0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩ 得出相应的点,t u （）； 4) 选择拉格朗日插值法，将,t u （）作为中间变量，在题目所给出的二维数表中进行二次代数插值，得到[][],)(z f x i y j =；5) 输出数表：[][][][]()()0,1,2,,10;0,1,2,,20,,,x i y j f x i y j i j =⋯=⋯； 6) 令k=0；7) 以()()(),,,0,1,r r r s x x y y r s ϕψ===…,k 为拟合基函数，将上述数表作为拟合条件，对于给定的k 值，得到矩阵B 、G 、U ；8) 令-1-1(),()T T T A B B B U C AG G G ==，用选主元的LU 分解法分别计算矩阵A 和C 的各列，最后得到系数矩阵C ；9) 以公式：()()()00,k ki j rs r i s j s r p x y C x y ϕψ===∑∑计算每个点的拟合值；10) 利用公式：()()()2102000,,i j i j i j f x y p x y σ===-∑∑计算拟合误差，当σ≤10-7时，循环结束，否则k=k+1，转（6）；11) 令[][]()**0.10.50.2 1,2,81,2,5x i i y j j i j ==+=⋯=⋯；，；，；12) 计算()()()******,,,,,i j i j i jf x y p x y delta x y ，输出数表，观察逼近效果； 13) 结束。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析大作业班级：学号：姓名：2015年11月5日注：给出最终报告的时候要有流程图，本人就是在报告中没有流程图所以大作业的分数并不高，最好讨论部分也要写的详细一些，之所以把FORTRAN语言版传上来是让用FORTRAN语言的同学在编程序时有个参照对比，还是建议大家自行编程，有助于加深对数值分析这门课的理解。
数值分析大作业第一题算法设计方案1.1分析s λ与1λ和501λ关系：题目中只给出了1λ和501λ的真实值大小关系，如果1λ和501λ同为正或者同为负，那么求出1λ和501λ后，在其中自然就找到s λ。
如果在求出1λ和501λ值后发现这两个特征值一正一负，则对原矩阵用一次反幂法即可得到s λ。
故可看出，本次任务的核心内容是要求出1λ和501λ的值。
1.2分析1λ和501λ的关系：假设1a λ=，501b λ=，先不考虑a b =这种情况① 0a b <≤，用幂法求出的值为a ，再用到原点平移原理，用幂法求出A aI -的按模最大特征值b a -，其中b a a -≤。
② 0a b ≤<，幂法求出值为b ，再用幂法求出A bI -的按摸最大特征值a b -，其中，a b b -<。
③ 0a b <<，当a b ≤时，幂法求出b ，再用幂法求出A bI -的按模最大特征值a b -，此时，a b b ->。
当a b >时，幂法求出a ，再用幂法求出A aI -的按模最大特征值b a -，此时b a a ->。
若a b =，则用幂法求出一值后，再求A 平移该值的按模最大特征值时会得到0这个结果。
1.3分析情况从而求解1λ和501λ：假设用幂法求出的值为x ，如果0x <，求A xI -的按模最大特征值y ，如果0y x <≤，对应情况①，此时1x λ=，501y x λ=+；如果0y =，则1501x λλ==，此时也可以认为1x λ=，501y x λ=+；如果y x >，此时1x λ=，501y x λ=+。
如果0x >，求A xI -的按模最大特征值y ，如果0y x <≤，对应情况②，此时501x λ=，1y x λ=+；如果0y =，则1501x λλ==，同上；如果y x >，501x λ=，1y x λ=+。
如果0x =，则15010λλ==。
求出1λ和501λ后，可根据1.1中的分析求出s λ。
2分析k μ从而求解ik λ：由于()501111,2,...,3940k kk λλμλ-=+=，在求出1λ和501λ后，即可判断k μ的正负。
在k 的变化过程中，可能k μ会从负变到正，也可能一直正负性一直不变。
无论k μ为正值还是负值，对k A I μ-用反幂法求得其按模最小值x ，ik k x λμ=+。
3求解()2cond A 和det A :由于()1222cond A AA -=，又由于A 是实对称矩阵，则T A A =，则2A ==2A 的全部特征值为正。
A 是实对称矩阵，则1A -也是实对称矩阵，又由于()221112A AAAA A I I ---===，则12A -===()1222cond A A A -==A 的任意特征值λ而言，Ax x λ=，()()22A x A x Ax x λλλ===，()()()max 2min A cond A A λλ=。
由于A 的最大特征值和最小特征值已在上面解出，易得()2cond A 。
()5011det det det det det ,I A LU L U U U I I =====∏。
注：在选择算法时，要求使矩阵A 的所有零元素都不储存，则采用一个二维数组()5,501C 存放A 中的带内元素即可，其中有转换关系21,ij i j j a c -++=。
源程序PROGRAM THE_FIRST_HOMEWORKIMPLICIT NONEINTEGER::M,N,R,S,IREAL(8)::X,Y,Z,BETA,ERROR,VALVE(1:3),CONDREAL(8),ALLOCATABLE::C(:,:),UNIT(:,:),MID(:,:),R_IK(:)OPEN(UNIT=12,FILE='DataOutput'//'\'//'1_The_Output_Of_The_First_Homework.txt') N=501M=5R=2S=2ERROR=1D-12ALLOCATE( C(1:501,1:5))ALLOCATE(UNIT(1:501,1:5))ALLOCATE( MID(1:501,1:5))ALLOCATE( R_IK(1:39))!把矩阵值放在二维数组中存储C(1:501,1)=-0.064D0;C(1 : 2,1)=0D0C(1:501,5)=-0.064D0;C(500:501,5)=0D0C(1:501,2)=0.16D0 ;C(1 : 1,2)=0D0C(1:501,4)=0.16D0 ;C(501:501,4)=0D0DO I=1,501C(I,3)=(1.64D0-0.024*DBLE(I))*DSIN(0.2*DBLE(I))-0.64*DEXP(0.1D0/DBLE(I))END DO!求R1,R501,|Rmin|CALL SOLVE_THE_EXTREME_FEATURE_VALVE(C,N,S,R,VALVE,ERROR)!求R_IKDO I=1,39MID=0D0R_IK(I)=VALVE(1)+I*(VALVE(2)-VALVE(1))/40D0MID(1:501,3)=R_IK(I)MID=C-MIDCALL MATRIX_FEATURE_VERSA_POWER_TWO_METHOD(MID,N,S,R,Z,ERROR)R_IK(I)=R_IK(I)+ZEND DO!求cond(A)2IF(DABS(VALVE(1))>DABS(VALVE(2)) )THENCOND=DABS(VALVE(1)/VALVE(3))ELSECOND=DABS(VALVE(2)/VALVE(3))END IF!先分解，之后求Det|A|CALL BREAK_DOWN_DOOLITTLE(C,N,S,R)X=1D0DO I=1,501X=X*C(I,3)END DOWRITE(12,*)'*************************计算结果*************************'WRITE(12,*)' 学号：姓名：WRITE(12,'("λ1 =",E19.12)')VALVE(1)WRITE(12,'("λ501 =",E19.12)')VALVE(2)WRITE(12,'("|λS| =",E19.12)')VALVE(3)WRITE(12,*)'λik ='WRITE(12,'(4(E19.12,2X))')R_IKWRITE(12,'("cond(A)2=",E19.12)')CONDWRITE(12,'("det|A| =",E19.12)')xWRITE(12,*)'*************************运算结束*************************'END PROGRAM THE_FIRST_HOMEWORK!求解一个矩阵的最大特征值和最小特征值和模最小特征值的程序，结果存储在VALVE中SUBROUTINE SOLVE_THE_EXTREME_FEATURE_VALVE(C,N,S,R,VALVE,ERROR)IMPLICIT NONEINTEGER::M,N,R,S,IREAL(8)::X,Y,ERROR,VALVE(1:3)REAL(8)::C(1:N,1:(S+R+1)),D(1:N,1:(S+R+1)),UNIT(1:N,1:(S+R+1))CALL MATRIX_FEATURE_POWER_INFINITE_METHOD(C,N,S,R,X,ERROR) !求C的按模最大特征值xUNIT=0D0UNIT(1:501,3)=XD=C-UNITCALL MATRIX_FEATURE_POWER_INFINITE_METHOD(D,N,S,R,Y,ERROR) !求C-XI的按模最大特征值yD=CIF(X<=0)THENVALVE(1)=X ;VALVE(2)=X+YIF(VALVE(2)<0)THENVALVE(3)=DABS(VALVE(2))ELSECALL MATRIX_FEATURE_VERSA_POWER_TWO_METHOD(D,N,S,R,VALVE(3),ERROR)VALVE(3)=DABS(VALVE(3))END IFELSEVALVE(1)=Y+X;VALVE(2)=XIF(VALVE(1)>0)THENVALVE(3)=VALVE(1)CALL MATRIX_FEATURE_VERSA_POWER_TWO_METHOD(D,N,S,R,VALVE(3),ERROR)VALVE(3)=DABS(VALVE(3))END IFEND IFRETURNEND SUBROUTINE SOLVE_THE_EXTREME_FEATURE_VALVE!求解矩阵的按模最大特征值和相应的特征向量，这里用的是幂法，用的是无穷范数的方法!C是二维数组，N是行长，M是列长，S是原矩阵的上半带宽，R是下半,BETA作为按模最大特!征值的存储，ERROR是误差限。
SUBROUTINE MATRIX_FEATURE_POWER_INFINITE_METHOD(C,N,S,R,BETA,ERROR)IMPLICIT NONEINTEGER N,M,I,J,K,REC,MAX_NUM,S,RREAL(8)::H,BETA0,BETA,ERRORREAL(8)::C(1:N,1:(S+R+1)),Y(1:N),U(1:N)!C存储原方程的各个系数，按列存储；Y存储逐次迭代的向量；U存储递推公式项!任意选U(0)作U的初步估计值,注意，当取U(1)=1D0,U(其它)=0D0时，结果不准确M=S+R+1U=1D0REC=0BETA=0D0MAX_NUM=1DOREC=REC+1BETA0=BETAH=DABS( U(MAX_NUM) )DO J=1,NIF(DABS( U(J) )>H)THENH=DABS( U(J) )MAX_NUM=JEND IFEND DOY=U/HDO J=1,NU(J)=0D0DO I=1,NIF((J-I+S+1)<1.OR.(J-I+S+1)>M)CYCLEU(J)=U(J)+C(I,J-I+S+1)*Y(I)END DOEND DOBETA=Y(MAX_NUM)*U(MAX_NUM)IF(REC>=2)THENIF( (DABS(BETA-BETA0)/DABS(BETA0))<=ERROR )EXITEND DORETURNEND SUBROUTINE MATRIX_FEATURE_POWER_INFINITE_METHOD!求解矩阵的按模最小特征值和对应的特征向量，这里用的是反幂法，并且用的是二范数的方法!C是二维数组，N是行长，M是列长，S是原矩阵的上半带宽，R是下半,BETA作为按模最小特!征值的存储，ERROR是误差限。