数值计算方法FORTRAN程序-函数插值和方程组求解

数值计算方法FORTRAN程序—函数插值和方程组求解吕毅宁吕毅宁：吕毅宁：lvyining@目录数值计算方法FORTRAN程序—函数插值和方程组求解 1一. 函数插值 1二. 齐次线性方程组的直接法求解 2三. 齐次线性方程组的迭代求解法 4附录一：流程图 7附录二：源程序 10程序一. Lagrange插值 10程序二. Newton向前插值 10程序三. Newton向后插值 11程序四. Gauss列主元消去法、Crout〔Dolittle〕三角分解法、LDLT、主程序 11程序五. Gauss列主元消去法子程序 12程序六. Doolittle三角分解子程序 12程序七. LDLT子程序 13程序八. 子程序 13程序九. Jacobi、Gauss-Seidel迭代法主程序 13程序十. Jacobi迭代法子程序 14程序十一. Gauss-Seidel迭代法子程序 14程序十二. SOR法主程序及子程序 15函数插值对n+1个节点xi及yi=f(xi)(I=0,1,…,n)，编制通用程序。n次Lagrange插值计算公式Ln(x);流程图〔图1、2〕程序〔程序一〕计算实例输入数据文件：input1.datKNOWNNODESNUMBERANDNODE'X'WITHY(X)WANTTOBEKNOWN:42.5NODESXIANDYIWHICHAREALREADYKNOWN:1.0003.0004.0002.0001.00027.00064.0008.000输出数据文件：output1.datResultsofLagrangeinterpolationatpointxis:Ln(x)=15.625000000000000结果分析实际计算数次插值，分析可见：在4节点插值时，对于四次以下的函数的插值具有很精确的结果，但对四次及以上的函数的插值，随阶次升高，插值精度降低。n次Newton向前插值计算公式;流程图〔图3〕程序〔程序二〕计算实例输入数据文件：input2.datKNOWNNODESNUMBERANDNODE'X'WITHY(X)WANTTOBEKNOWN:42.5NODESXIANDYIWHICHAREALREADYKNOWN:1.0003.0004.0002.0001.00081.000256.00016.000FLAGTOCONTROLWHETHERTOPRINTDIFF-CHART:1输出数据文件：output2.datNn(x)=38.500000000000000Forwarddifferencediagram:1.00000000016.0000000015.0000000081.0000000065.0000000050.00000000256.0000000175.0000000110.000000060.00000000结果分析本算例用四个节点的Newton向前插值计算公式对四次函数进行插值，计算误差为：绝对误差： |38.500000000000000-2.54|=|38.5-39.0625|=0.5625;相对误差: 0.5625/39.0625=1.44%;n次Newton向后插值计算公式;流程图〔与n次Newton向前插值计算时类似〕程序〔程序三〕计算实例输入数据文件：input3.datKNOWNNODESNUMBERANDNODE'X'WITHY(X)WANTTOBEKNOWN:42.5NODESXIANDYIWHICHAREALREADYKNOWN:1.0003.0004.0002.0001.00027.00064.0008.000FLAGTOCONTROLWHETHERTOPRINTDIFF-CHART:1输出数据文件：output3.datNn(x)=15.625000000000000Backwarddifferencediagram:64.0000000027.00000000-37.000000008.000000000-19.0000000018.000000001.000000000-7.00000000012.00000000-6.000000000结果分析本算例用四个节点的Newton向后插值计算公式对三次函数进行插值，计算结果精确。齐次线性方程组的直接法求解编制解Ax=b的通用子程序。列主元消去法；流程图〔图4、5〕程序〔程序四、五〕计算实例输入数据文件input.dat：COEFFICIENTSMATRIXIS:100.0.0.37.4.0.46.400.0.110.0.38.0.0.200.0.40.55.73.0.0.100.61.0.20.19.-0.-0.100.0.0.75.0.0.0.200.RIGHTVECTORIS:20.11.23.63.72.91.输出数据文件output.datAnswertothesimutaneouslinearequationsis:-4.312691791318356E-001-6.207237146208434E-0021.150000000000000E-0018.578492048870718E-0011.679627521841597E-0014.351687505777960E-001VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:3.552713678800501E-0157.105427357601002E-0150.000000000000000E+0001.421085471520200E-0141.065814103640150E-0140.000000000000000E+000结果分析从output.dat中的校验数据〔VERIFYDATA〕可见，方程组的求解精度可以到达1.0*10-14〔校验数据=原来右端向量-原来系数矩阵*求解结果〕。实现矩阵三角分解的Doolittle及Cholesky三角分解方法及用此方法解Ax=b的过程；Doolittle三角分解方法流程图〔图6〕程序〔程序四、六〕计算实例输入数据文件input.dat：〔与上面列主元消去法时相同〕输出数据文件output.dat：Answertothesimutaneouslinearequationsis:-4.312691791318356E-001-6.207237146208434E-0021.150000000000000E-0018.578492048870718E-0011.679627521841597E-0014.351687505777960E-001VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:3.552713678800501E-0157.105427357601002E-0150.000000000000000E+0001.421085471520200E-0141.065814103640150E-0140.000000000000000E+000结果分析从output.dat中的校验数据〔VERIFYDATA〕可见，方程组的求解精度可以到达1.0*10-15.与上面列主元消去法时的计算结果比拟，可见有相同的很高的求解精度。Cholesky三角分解方法〔见下面4〕〕用列主元消去法程序解方程组并比拟计算结果精度〔准确解为x1=x2=x3=x4=1〕输入数据文件input.dat:COEFFICIENTSMATRIXIS:1.13480.53013.41291.23713.83261.78574.93174.99981.16512.53308.764310.67213.40171.54351.31420.0147RIGHTVECTORIS:9.53246.394118.423116.9237输出数据文件output.dat:Answertothesimutaneouslinearequationsis:9.997246338638249E-0019.971273542063148E-0011.0013745262416351.002328461426957VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:1.332267629550188E-0156.661338147750939E-0161.776356839400251E-0152.726985304235541E-015结果分析从output.dat中的校验数据〔VERIFYDATA〕可见，方程组的求解精度可以到达1.0*10-15.用平方根法解线性方程组流程图〔与Doolittle三角分解类似〕程序〔程序四、七〔LDLT〕、八〔〕〕计算实例输入数据文件input.datCOEFFICIENTSMATRIXIS:4.2.42.3.2.45.444.5.82.4.5.217.453.5.87.4519.66RIGHTVECTORIS:12.28016.92822.95750.945输出数据文件output.datAnswertothesimutaneouslinearequationsis:1.200000000000000-7.999999999999994E-0011.7000000000000002.000000000000000VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:8.881784197001252E-0160.000000000000000E+0003.552713678800501E-0150.000000000000000E+000结果分析从output.dat中的校验数据〔VERIFYDATA〕可见，方程组的求解精度可以到达1.0*10-15.齐次线性方程组的迭代求解法编制解Ax=b的通用子程序。Jacobi

迭代格式；流程图〔图7〕程序〔程序十〕计算实例输入数据文件input.datMAXIMUMITERATIONNUMBERANDEPSISONARE:100001.0D-6COEFFICIENTSMATRIXIS:1.00.-.25-.250.1.00-.25-.25-.25-.251.000.-.25-.250.1.00RIGHTVECTORIS:0.50.50.50.5输出数据文件output.datAnswertothesimultaneouslinearequationsis:9.999990463256836E-0019.999990463256836E-0019.999990463256836E-0019.999990463256836E-001VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:4.768371582031250E-0074.768371582031250E-0074.768371582031250E-0074.768371582031250E-007迭代次数：20结果分析对课本上的例题实际用Jacobi迭代格式计算，迭代20次后收敛到满足给定精度的近似解。Gauss-Seidel迭代格式；流程图〔与Jacobi迭代格式时类似〕程序〔程序十一〕计算实例输入数据文件input.dat〔与上面Jacobi

迭代格式时相同〕输出数据文件output.datAnswertothesimultaneouslinearequationsis:9.999997615814209E-0019.999997615814209E-0019.999998807907105E-0019.999998807907105E-001VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:1.788139343261719E-0071.788139343261719E-0070.000000000000000E+0000.000000000000000E+000迭代次数：11结果分析对上面的算例用Gauss-Seidel迭代格式计算，迭代11次后收敛到满足给定精度的近似解。与用Jacobi迭代格式计算时相比，收敛速度更快。SOR方法格式流程图〔与Jacobi迭代格式时类似〕程序〔程序十二〕计算实例输入数据文件input.datMAXIMUMITERATIONNUMBERANDEPSISONARE:100001.0D-6 1.1D+0(方程组系数及右端向量与Jacobi迭代格式时相同。)输出数据文件output.datAnswertothesimutaneouslinearequationsis:9.999995649140370E-0019.999995649140370E-0019.999998505886739E-0019.999998505886739E-001VERIFYDATA.DifferencebetweentheOriginalVectorandMatrix*P:3.603802999907479E-0073.603802999907479E-0076.813165531749377E-0086.813165531749377E-008迭代次数：7结果分析令松弛因子等距变化，对相同精度要求，迭代次数如下：松弛因子0.91.01.41.5迭代次数1612710131722对此例，取=1.1时所用迭代次数最少。可见，适中选取可得到比Jacobi迭代格式和Gauss-Seidel迭代格式更快的收敛速度。附录一：流程图〔图1〕n次Lagrange插值-流程图 〔图2〕Lagrange插值模块流程图参数说明：入口：X： 被插值点；XI，YI： 插值点的坐标；N： 插值多项式次数；出口： XLNX： X点的插值函数值；SUBLAGINTPO子程序流程：XLNX=0.0D0，K=0；I=0，XLKX=1.0D0； y I=?Kn XLKX=XLKX*(X-XI(I))/ (XI(K)-XI(I)) I>N XLNX=XLNX+XLKX*YI(K)K>NXLNX=0.0D0，K=0；I=0，XLKX=1.0D0； y I=?Kn XLKX=XLKX*(X-XI(I))/ (XI(K)-XI(I)) I>N XLNX=XLNX+XLKX*YI(K)K>NJ=J+1NI=I+1NI=2YI>N+1YYNJ>II>N+1子程序结束TMP=TMP*(T-I+2)/(I-1)IC=IC+IXNNX=XNNX+TMP*XTEMP(IC)T=N*(X-XI(0))/(XI(N)-XI(0))IC=1TMP=1.0D0XNNX=XTEMP(IC)I=1XTEMP(IC)=XTEMP(IC-1)-XTEMP(IC-I)IC=IC+1J=2IC=1XTEMP(IC)=YI(I-1)IC=IC+1 〔图3〕Newton向前插值子程序流程图J=J+1NI=I+1NI=2YI>N+1YYNJ>II>N+1子程序结束TMP=TMP*(T-I+2)/(I-1)IC=IC+IXNNX=XNNX+TMP*XTEMP(IC)T=N*(X-XI(0))/(XI(N)-XI(0))IC=1TMP=1.0D0XNNX=XTEMP(IC)I=1XTEMP(IC)=XTEMP(IC-1)-XTEMP(IC-I)IC=IC+1J=2IC=1XTEMP(IC)=YI(I-1)IC=IC+1程序终止插值结果输出SUBOUTPUTLagrange插值计算SUBLAGINTPO数据升序排列SUBSORTXI输入数据SUBINPUT开始程序终止插值结果输出SUBOUTPUTLagrange插值计算SUBLAGINTPO数据升序排列SUBSORTXI输入数据SUBINPUT开始I=I+1

〔图4〕列主元Gauss消去法程序流程图 〔图5〕Gauss求解子模块流程图I=I+1程序终止结果校验SUBVERIFY求解结果输出SUBOUTPUT求解方程组SUBGSSOLVER备份系数矩阵A和右端项P输入数据SUBINPUT开始程序终止结果校验SUBVERIFY求解结果输出SUBOUTPUT求解方程组SUBGSSOLVER备份系数矩阵A和右端项P输入数据SUBINPUT开始参数说明：入口：A： 系数矩阵；P： 方程组右端项；N： 方程组阶数；出口： P： 方程组的解；SUBGSSOLVER子程序流程：寻找列主元，CALLDOMCOLNYI=NYK>NJ=J+1A(I,J)=A(I,J)-A(I,K)*A(K,J)J>NNYSTOPN|A(K,K)|<1.0E-6I=K+1NNJ=J+1NYJ>NK=1子程序结束I<1P(I)=P(I)/A(I,I)P(I)=P(I)-A(I,J)*P(J)J=I+1J=K+1A(I,K)=A(I,K)/A(K,K)P(I)=P(I)-A(I,K)*P(K)寻找列主元，CALLDOMCOLNYI=NYK>NJ=J+1A(I,J)=A(I,J)-A(I,K)*A(K,J)J>NNYSTOPN|A(K,K)|<1.0E-6I=K+1NNJ=J+1NYJ>NK=1子程序结束I<1P(I)=P(I)/A(I,I)P(I)=P(I)-A(I,J)*P(J)J=I+1J=K+1A(I,K)=A(I,K)/A(K,K)P(I)=P(I)-A(I,K)*P(K)Gauss消去法程序说明：为了程序的简洁，易于调试，没有开大数组进行动态内存分配，来提高主程序的通用性；Gauss消去法程序说明：为了程序的简洁，易于调试，没有开大数组进行动态内存分配，来提高主程序的通用性；为求解不同规模的方程组，需要在主程序中修改参数；为了验证求解结果，另开数组备份了系数矩阵和右端项；K=K+1

〔图6〕三角分解Doolittle法子模块流程图NYNNYI=I+1NNJ=J+1NI=I+1J=J+1J>I-1I>NYYI=NJ=I+1P(I)=P(I)-A(I,J)*P(J)J>NP(I)=P(I)-A(I,J)*P(J)I<1Y子程序结束NK=K+1YI>NYI=I+1J=J+1A(K,I)=A(K,I)/A(I,I)YI=1YP(I)=P(I)-A(I,J)*P(J)K>NI=1J=1J>K-1A(K,I)=A(K,I)-A(K,J)*A(J,I)J=1J>I-1I>K-1A(K,I)=A(K,I)-A(K,J)*A(J,I)J=1I=1K=1NYNNYI=I+1NNJ=J+1NI=I+1J=J+1J>I-1I>NYYI=NJ=I+1P(I)=P(I)-A(I,J)*P(J)J>NP(I)=P(I)-A(I,J)*P(J)I<1Y子程序结束NK=K+1YI>NYI=I+1J=J+1A(K,I)=A(K,I)/A(I,I)YI=1YP(I)=P(I)-A(I,J)*P(J)K>NI=1J=1J>K-1A(K,I)=A(K,I)-A(K,J)*A(J,I)J=1J>I-1I>K-1A(K,I)=A(K,I)-A(K,J)*A(J,I)J=1I=1K=1Y子程序结束为方程组的解判断是否收敛CALLCONVERGNITRA=ITRA+1ITRA=0提供迭代初值，P(I)=0.0，(I=1,N)由BJ=D-1(L+U),及fJ=D-1b,求BJ和fJY子程序结束为方程组的解判断是否收敛CALLCONVERGNITRA=ITRA+1ITRA=0提供迭代初值，P(I)=0.0，(I=1,N)由BJ=D-1(L+U),及fJ=D-1b,求BJ和fJ〔图7〕Jocabi迭代法子模块流程图附录二：源程序Lagrange插值 PROGRAMLAGRANGE_INTERPOLATION IMPLICITREAL*8(A-H,O-Z) PARAMETER(NMAX=3000) PARAMETER(IU1=8,IU2=9) DIMENSIONXI(0:NMAX),YI(0:NMAX) CALLINPUT(IU1,XI,YI,X,N,NMAX) CALLSORTXI(XI,YI,N) CALLLAGINTPO(XLNX,X,XI,YI,N) CALLOUTPUT(IU2,XLNX) STOP'Exitwithcode(0)...' END SUBROUTINEINPUT(IU1,XI,YI,X,N,NMAX) IMPLICITREAL*8(A-H,O-Z) DIMENSIONXI(0:NMAX),YI(0:NMAX) CHARACTER*80AAAA WRITE(*,*)'PleaseinputanINPUTfilename:' READ(*,*)AAAA OPEN(IU1,FILE=AAAA,STATUS='UNKNOWN') READ(IU1,'(80A)')AAAA READ(IU1,*)N,X N=N-1 READ(IU1,'(80A)')AAA READ(IU1,*)(XI(I),I=0,N) READ(IU1,*)(YI(I),I=0,N) RETURN END SUBROUTINEOUTPUT(IU2,XLNX) IMPLICITREAL*8(A-H,O-Z) CHARACTER*80AAAA WRITE(*,*)'PleaseinputanOUTPUTfilename:' READ(*,*)AAAAOPEN(IU2,FILE=AAAA,status='UNKNOWN') WRITE(IU2,*)'ResultsofLagrangeinterpolationatpointxis:' WRITE(IU2,*)'Ln(x)=',XLNXRETURN END SUBROUTINELAGINTPO(XLNX,X,XI,YI,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONXI(0:N),YI(0:N) XLNX=0.0D0 DO10,K=0,N CALLPRODUCT(XLKX,X,K,XI,N) XLNX=XLNX+XLKX*YI(K)10 CONTINUERETURN END SUBROUTINEPRODUCT(XLKX,X,K,XI,N) IMPLICITREAL*8(A-H,O-Z)DIMENSIONXI(0:N) XLKX=1.0D0 DO20,I=0,N IF(I.EQ.K)GOTO20 XLKX=XLKX*(X-XI(I))/(XI(K)-XI(I))20 CONTINUERETURN END SUBROUTINESORTXI(XI,YI,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONXI(0:N),YI(0:N) DO10,I=0,N-1 DO20,J=I+1,N IF(XI(J).LT.XI(I))THEN TMP=XI(J) XI(J)=XI(I) XI(I)=TMP TMP=YI(J) YI(J)=YI(I) YI(I)=TMP ENDIF20 CONTINUE10 CONTINUERETURN ENDNewton向前插值 PROGRAMNEWTON_FORWARD_INTERPOLATION IMPLICITREAL*8(A-H,O-Z) PARAMETER(NM=1000) PARAMETER(IU1=8,IU2=9) DIMENSIONXI(0:NM),YI(0:NM),XTEMP((NM+1)*(NM+2)/2) CALLINPUT(IU1,XI,YI,X,N,NPC,NM) CALLSORTXI(XI,YI,N) CALLNTFINTPO(XNNX,X,XI,YI,XTEMP,N) CALLOUTPUT(IU2,XNNX,XTEMP,N,NPC) STOP'Exitwithcode(0)...' END SUBROUTINEINPUT(IU1,XI,YI,X,N,NPC) IMPLICITREAL*8(A-H,O-Z)DIMENSIONXI(0:NM),YI(0:NM) OPEN(IU1,FILE='Input.dat.ntf.4',STATUS='OLD') READ(IU1,'(80A)')AAA READ(IU1,*)N,X N=N-1 READ(IU1,'(80A)')AAA READ(IU1,*)(XI(I),I=0,N) READ(IU1,*)(YI(I),I=0,N) READ(IU1,'(80A)')AAA READ(IU1,*)NPC RETURN END SUBROUTINEOUTPUT(IU2,XNNX,XTEMP,N,NPC) IMPLICITREAL*8(A-H,O-Z)DIMENSIONXTEMP((N+1)*(N+2)/2) OPEN(IU2,FILE='Output.dat.ntf',STATUS='UNKNOWN') WRITE(IU2,*)'Nn(x)=',XNNX IF(NPC.NE.0)THEN WRITE(IU2,*)'Forwarddifferencediagram:' IC=0 DO10,I=1,N+1 WRITE(IU2,'(<N+1>G20.10)')(XTEMP(J),J=IC+1,IC+I) IC=IC+I10CONTINUEENDIFRETURN ENDSUBROUTINENTFINTPO(XNNX,X,XI,YI,XTEMP,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONXI(0:N),YI(0:N),XTEMP((N+1)*(N+2)/2)IC=1 DO10,I=1,N+1 XTEMP(IC)=YI(I-1) IC=IC+1 DO20,J=2,I XTEMP(IC)=XTEMP(IC-1)-XTEMP(IC-I) IC=IC+120 CONTINUE10 CONTINUET=N*(X-XI(0))/(XI(N)-XI(0)) IC=1 TMP=1.0D0 XNNX=XTEMP(IC) DO30,I=2,N+1 TMP=TMP*(T-I+2)/(I-1) IC=IC+I XNNX=XNNX+TMP*XTEMP(IC)30 CONTINUERETURN END SUBROUTINESORTXI(XI,YI,N)〔同程序一〕Newton向后插值 PROGRAMNEWTON_BACKWARD_INTERPOLATION IMPLICITREAL*8(A-H,O-Z) PARAMETER(N=3) PARAMETER(IU1=8,IU2=9) DIMENSIONXI(0:N),YI(0:N),XTEMP((N+1)*(N+2)/2) CALLINPUT(IU1,XI,YI,X,N,NPC) CALLSORTXI(XI,YI,N) CALLNTBINTPO(XNNX,X,XI,YI,XTEMP,N) CALLOUTPUT(IU2,XNNX,XTEMP,N,NPC) STOP'Exitwithcode(0)...' END SUBROUTINEINPUT(IU1,XI,YI,X,N,NPC)〔同程序二〕 SUBROUTINEOUTPUT(IU2,XNNX,XTEMP,N,NPC)〔同程序二〕SUBROUTINENTBINTPO(XNNX,X,XI,YI,XTEMP,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONXI(0:N),YI(0:N),XTEMP((N+1)*(N+2)/2)IC=1 DO10,I=1,N+1 XTEMP(IC)=YI(N+1-I) IC=IC+1 DO20,J=2,I XTEMP(IC)=XTEMP(IC-1)-XTEMP(IC-I) IC=IC+120 CONTINUE10 CONTINUET=N*(X-XI(0))/(XI(N)-XI(0)) IC=1 TMP=1.0D0 XNNX=XTEMP(IC) DO30,I=2,N+1 TMP=TMP*(T-I+2)/(I-1) IC=IC+I XNNX=XNNX+TMP*XTEMP(IC)30 CONTINUERETURN END SUBROUTINESORTXI(XI,YI,N)〔同程序二〕Gauss列主元消去法、Crout〔Dolittle〕三角分解法、LDLT、主程序 PROGRAMGAUSS_SOLVER IMPLICITREAL*8(A-H,O-Z)PARAMETER(N=6) PARAMETER(IU1=8,IU2=9) DIMENSIONA(N,N),P(N) DIMENSIONB(N,N),Q(N) CALLINPUT(IU1,A,P,N)DO10,I=1,N Q(I)=P(I) DO10,J=1,N B(J,I)=A(J,I)10CONTINUE CALLGSSOLVER(A,P,N) CALLOUTPUT(IU2,P,N)CALLVERIFYI(IU2,P,Q,B,N) STOP'Exitwithcode(0)...' END SUBROUTINEVERIFYI(IU2,P,Q,B,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONP(N),B(N,N),Q(N) DO10,I=1,N DO10,K=1,N Q(I)=Q(I)-B(I,K)*P(K)10 CONTINUEWRITE(IU2,*)'VERIFYDATA.' WRITE(IU2,100)WRITE(IU2,*)(DABS(Q(I)),I=1,N)100 FORMAT('DifferencebetweentheOriginalVectorandMatrix*P:') RETURN END SUBROUTINEINPUT(IU1,A,P,N) IMPLICITREAL*8(A-H,O-Z) CHARACTER*80AAA DIMENSIONA(N,N),P(N)OPEN(IU1,FILE='LynSolverInput1.dat',STATUS='OLD') READ(IU1,*)AAA READ(IU1,*)A READ(IU1,*)AAAREAD(IU1,*)P RETURN END SUBROUTINEOUTPUT(IU2,P,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONP(N)OPEN(IU2,FILE='LynSolverOutput.dat') WRITE(IU2,*)'Answertothesimultaneouslinearequationsis:'WRITE(IU2,*)(P(I),I=1,N) RETURN ENDGauss列主元消去法子程序SUBROUTINEGSSOLVER(A,P,N) IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N) DO10,K=1,N !column CALLDOMCOL(A,P,N,K)IF(DABS(A(K,K)).LT.1.0D-6)STOP'Currentmatrixissingular...' DO10,I=K+1,N !row A(I,K)=A(I,K)/A(K,K) P(I)=P(I)-A(I,K)*P(K) DO30,J=K+1,N A(I,J)=A(I,J)-A(I,K)*A(K,J)30 CONTINUE10 CONTINUE DO20,I=N,1,-1 DO40,J=I+1,N P(I)=P(I)-A(I,J)*P(J)40 CONTINUE P(I)=P(I)/A(I,I)20 CONTINUERETURN END SUBROUTINEDOMCOL(A,P,N,K) IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N) TEMP=DABS(A(K,K)) ITEMP=K DO10,I=K+1,N IF(DABS(A(I,K)).GT.TEMP)THEN TEMP=A(I,K) ITEMP=I ENDIF10CONTINUEIF(ITEMP.NE.K)THEN DO20,I=K,N TEMP=A(K,I) A(K,I)=A(ITEMP,I) A(ITEMP,I)=TEMP20 CONTINUE TEMP=P(K) P(K)=P(ITEMP) P(ITEMP)=TEMPENDIFRETURN ENDDoolittle三角分解子程序SUBROUTINEDOOLITTLESOLVER1(A,P,N)!ROWBYROW; IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N) DOK=1,N DOI=1,K-1 DOJ=1,I-1 A(K,I)=A(K,I)-A(K,J)*A(J,I) ENDDO A(K,I)=A(K,I)/A(I,I) ENDDO DOI=K,N DOJ=1,K-1 A(K,I)=A(K,I)-A(K,J)*A(J,I) ENDDO ENDDO ENDDO DOI=1,N DOJ=1,I-1 P(I)=P(I)-A(I,J)*P(J) ENDDO ENDDO DOI=N,1,-1 DOJ=I+1,N P(I)=P(I)-A(I,J)*P(J) ENDDO P(I)=P(I)/A(I,I) ENDDO RETURN ENDSUBROUTINEDOOLITTLESOLVER2(A,P,N)!COLUMNBYCOLUMN; IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N) DOJ=1,N DOI=1,J DOK=1,I-1 A(I,J)=A(I,J)-A(I,K)*A(K,J) ENDDO ENDDO DOI=J+1,N DOK=1,J-1 A(I,J)=A(I,J)-A(I,K)*A(K,J) ENDDO A(I,J)=A(I,J)/A(J,J) ENDDO ENDDO DOI=1,N DOJ=1,I-1 P(I)=P(I)-A(I,J)*P(J) ENDDO ENDDO DOI=N,1,-1 DOJ=I+1,N P(I)=P(I)-A(I,J)*P(J) ENDDO P(I)=P(I)/A(I,I) ENDDO RETURN ENDLDLT子程序SUBROUTINELDLTSOLVER(A,P,N) !BOOK.p.179 IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N) DOI=1,N DOJ=1,I IF(DABS(A(I,J)-A(J,I)).GT.1.0D-3)WRITE(*,*)'*Aisnotsymmetric*' ENDDO ENDDO DOK=1,N DOJ=1,K-1 DOI=1,J-1 A(J,K)=A(J,K)-A(J,I)*A(I,K) ENDDO A(J,K)=A(J,K)/A(J,J) ENDDO DOJ=K,N DOI=1,K-1 A(J,K)=A(J,K)-A(J,I)*A(I,K) ENDDO ENDDO ENDDO DOI=1,N DOJ=I+1,N A(J,I)=A(J,I)/A(I,I) ENDDO ENDDO DOI=1,N DOJ=1,I-1 P(I)=P(I)-A(I,J)*P(J) ENDDO ENDDO DOI=N,1,-1 P(I)=P(I)/A(I,I) DOJ=I+1,N P(I)=P(I)-A(I,J)*P(J) ENDDO ENDDO RETURN END子程序SUBROUTINELLTSOLVER(A,P,N)!BOOK.p.176 IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N) DOJ=1,N DOK=1,J-1 A(J,J)=A(J,J)-A(J,K)*A(K,J) ENDDO A(J,J)=DSQRT(A(J,J)) DOI=J+1,N DOK=1,J-1 A(I,J)=A(I,J)-A(I,K)*A(K,J) ENDDO A(I,J)=A(I,J)/A(J,J) A(J,I)=A(I,J) ENDDO ENDDO DOI=1,N DOK=1,I-1 P(I)=P(I)-A(I,K)*P(K) ENDDO P(I)=P(I)/A(I,I) ENDDO DOI=N,1,-1 DOK=I+1,N P(I)=P(I)-A(I,K)*P(K) ENDDO P(I)=P(I)/A(I,I) ENDDO RETURN ENDJacobi、Gauss-Seidel迭代法主程序 PROGRAMJACOB_SOLVER IMPLICITREAL*8(A-H,O-Z) PARAMETER(N=4) DIMENSIONA(N,N),P(N),P0(N),PP(N) DIMENSIONB(N,N),Q(N) CALLINPUT(IU1,A,P,N,EPS,MITRA)DO10,I=1,N Q(I)=P(I) DO10,J=1,N B(J,I)=A(J,I)10CONTINUECALLJACOBSOLVER(A,P,P0,PP,N,EPS,MITRA) CALLOUTPUT(IU2,P,N)CALLVERIFYI(IU2,P,Q,B,N) STOP'Exitwithcode(0)...' END SUBROUTINEVERIFYI(IU2,P,Q,B,N)〔同程序三〕 SUBROUTINEINPUT(IU1,A,P,N,EPS,MITRA) IMPLICITREAL*8(A-H,O-Z) CHARACTER*80AAA DIMENSIONA(N,N),P(N)OPEN(IU1,FILE='LynJacobSolverInput2.dat',STATUS='OLD') READ(IU1,*)AAA READ(IU1,*)MITRA,EPS READ(IU1,*)AAA READ(IU1,*)A READ(IU1,*)AAAREAD(IU1,*)P RETURN END SUBROUTINEOUTPUT(IU2,P,N)〔同程序四〕Jacobi迭代法子程序〔主程序与程序四类似〕SUBROUTINEJACOBSOLVER(A,P,P0,PP,N,EPS,MITRA) IMPLICITREAL*8(A-H,O-Z) DIMENSIONA(N,N),P(N),P0(N),PP(N) DOI=1,N DOJ=1,N A(J,I)=-A(J,I) ENDDO A(I,I)=-A(I,I) ENDDO DOI=1,N PP(I)=P(I)/A(I,I) DOJ=1,I-1 A(I,J)=A(I,J)/A(I,I) ENDDO DOJ=I+1,N A(I,J)=A(I,J)/A(I,I) ENDDO ENDDO DOI=1,N A(I,I)=0.0D0 ENDDO DOI=1,N P0(I)=0.0D0 ENDDO ITRA=0100 CONTINUE DOI=1,N P(I)=0.0D0 DOJ=1,N P(I)=P(I)+A(I,J)*P0(J) ENDDO P(I)=P(I)+PP(I) ENDDO IFLAG=0 CALLCONVERG(IFLAG,ITRA,P,P0,N,EPS,MITRA) DOI=1,N P0(I)=P(I) ENDDO IF(IFLAG.EQ.0)GOTO100 RETURN END SUBROUTINECONVERG(IFLAG,ITRA,P,P0,N,EPS,MITRA) IMPLICITREAL*8(A-H,O-Z) DIMENSIONP(N),P0(N) ITRA=ITRA+1 write(*,*)'itera