正交配置求解问题.doc

正交配置求解问题 正交配置求解问题 运用正交配置法求解有轴向扩散的固定床反应器中催化反应的温度和浓度分布 柱形固体 床反应器中催化反应的温度和浓度方程为 Z T r 1 r T r r TcR Z c r c r rr 1 TcR r 1 r 0 0 r T r c r 1 Biw T 1 z Tw z r 1 0 r T r c T r 0 T0 c r 0 c0 Z T Z T r 1 r T r r TcR Z c r c r rr 1 TcR r 1 r 0 0 r T r c r 1 Biw T 1 z Tw z r 1 0 r T r c T r 0 T0 c r 0 c0 其中 R c T 为催化反应的速率方程 其形式为 R c T T ec 1 r 1 r T r r TcR Z c r c r rr 1 TcR r 1 r 0 0 r T r c r 1 Biw T 1 z Tw z r 1 0 r T r c T r 0 T0 c r 0 c0 其中 R c T 为催化反应的速率方程 其形式为 R c T T ec 1 解题思路 解题思路 应用对称的正交配置法 有下面的方程和初始条件 1 dZ j dT i T ji B N i 1 1 j c T e 1 dZ j dc i c ji B N i 1 1 j c T e Tj 0 T0 cj 0 c0 边界条件为 AN 1 iTi Biw TN 1 Tw AN 1 ici 0 1 1 N i 1 1 N i 将温度和浓度的边界条件代入微分方程 消去边界值 可得 2N 个常微分方程 而将两边界 条件的代数方程同 2N 个常微分方程组联合 就组成 2N 2 个微分代数方程组 结合正交配 置系数的计算程序与常微分方程组或微分方程组求解程序 可得到反应器中的温度和浓度 分布 具体做法如下 一 利用对称的正交配置格式 1 对称常微分方程程序 COLLAB FOR DLSODE FOR 主程序 IMPLICIT REAL 8 A H O Z EXTERNAL FEX JEX DIMENSION AS 19 19 BS 19 19 Q 19 19 XS 19 WS 19 DIMENSION DIF1 19 DIF2 19 DIF3 19 ROOT 19 V1 19 V2 19 DIMENSION Y 99 ATOL 99 RWORK 10920 IWORK 120 DOUBLE PRECISION YN1 YN2 COMMON AB N AS BS COMMON BC YN1 YN2 C N FOR SYMMETRIC COLLOCATION USED FOR PARTICLE AND C M FOR ASYMMETRIC COLLOCATION USED FOR COLUMN N 7 IW 1 IS 2 CALL COLL AS BS Q XS WS 19 N IW IS NS N 1 WRITE SYMMETRIC SITUATION WRITE POLYNOMIAL ROOTS WRITE XS I I 1 NS WRITE WRITE A MATRIX DO 20 I 1 NS 20 WRITE AS I J J 1 NS WRITE WRITE B MATRIX DO 30 I 1 NS 30 WRITE BS I J J 1 NS WRITE WRITE W MATRIX WRITE WS J J 1 NS C CALCULATING THE PARAMETERS OF THE PROBLEM WHICH WILL BE USED C FOR THE DIMENSIONLESS FORM OF AND DEFINING OF THE PROBLEM NEQ 2 N LRW 22 9 NEQ NEQ 2 LIW 20 NEQ C INITIAL CONDITIONS DO 201 I 1 N Y I 1 D0 Y N I 0 D0 201 CONTINUE YN1 1 0D0 YN2 0 D0 T 0 D0 DT 5 D 2 ITOL 2 RTOL 1 D 6 DO 203 I 1 NEQ ATOL I 1 D 6 203 CONTINUE ITASK 1 ISTATE 1 IOPT 0 MF 22 DO 240 IOUT 1 20 TOUT DT DFLOAT IOUT CALL LSODE FEX NEQ Y T TOUT ITOL RTOL ATOL ITASK ISTATE 1 IOPT RWORK LRW IWORK LIW JEX MF OPEN 2 FILE LW S ODE OUT WRITE 2 Z F8 4 T WRITE 2 R WRITE 2 10 4X D11 5 XS I I 1 N 1 WRITE 2 T WRITE 2 10 4X D11 5 Y I I 1 N YN1 WRITE 2 C C DO 205 I 1 N WRITE 2 10 4X D11 5 Y N I I 1 N YN2 C WRITE 2 C 205 CONTINUE C 220 FORMAT 7H AT T D12 4 6H Y 3D15 7 IF ISTATE LT 0 GO TO 280 240 CONTINUE WRITE 3 260 IWORK 11 IWORK 12 IWORK 13 260 FORMAT 12H NO STEPS I4 11H NO F S I4 11H NO J S I4 STOP 280 WRITE 3 290 ISTATE 290 FORMAT 22H ERROR HALT ISTATE I3 STOP END 子程序 SUBROUTINE FEX NEQ T Y YP C VARIABLES IN THE ORDER OF Y 1 TO Y M FOR THE COLUMN COLLOCATION C POINTS OF 2 TO M 1 Y M 1 TO Y M N FOR PARTICLE COLLOCATION POINTS C FROM 1 TO N IN COLUMN COLLOCATION POINT 1 AND Y M N M 1 1 TO C Y M N M 2 FOR PARTICLE COLLOCATION POINTS FROM 1 TO N IN COLUMN C COLLOCATION POINT M 2 IMPLICIT REAL 8 A H O Z DIMENSION Y 19 19 YP 19 19 AS 19 19 BS 19 19 DOUBLE PRECISION YN1 YN2 COMMON AB N AS BS COMMON BC YN1 YN2 T0 0 D0 TT0 0 D0 DO 990 I 1 N T0 T0 AS N 1 I Y I TT0 TT0 AS N 1 I Y N I 990CONTINUE YN1 0 92 T0 1 AS N 1 N 1 YN2 TT0 AS N 1 N 1 DO 99 J 1 N T1 BS J N 1 YN1 T2 BS J N 1 YN2 DO 991 I 1 N T1 T1 BS J I Y I T2 T2 BS J I Y I N 991CONTINUE YP J T1 0 2 1 Y N J EXP 20 20 Y J YP N J T2 0 3 1 Y N J EXP 20 20 Y J 99 CONTINUE RETURN END SUBROUTINE JEX NEQ T Y ML MU PD NRPD DOUBLE PRECISION PD T Y DIMENSION Y NEQ PD NRPD NEQ RETURN END 输出结果 Z 0500 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 10109D 01 10104D 01 10088D 01 10056D 01 10008D 01 99565D 00 99154D 00 98955D 00 C 16524D 01 16408D 01 16087D 01 15522D 01 14890D 01 14453D 01 14290D 01 14878D 01 Z 1000 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 10215D 01 10192D 01 10150D 01 10091D 01 10025D 01 99642D 00 99205D 00 99015D 00 C 35719D 01 34822D 01 33344D 01 31668D 01 30296D 01 29520D 01 Z 9500 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 13623D 01 13464D 01 13210D 01 12905D 01 12599D 01 12338D 01 12155D 01 12046D 01 C 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 Z 1 0000 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 13290D 01 13142D 01 12906D 01 12624D 01 12341D 01 12099D 01 11930D 01 11825D 01 C 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 NO STEPS 210 NO F S 677 NO J S 29 2 对称代数微分方程程序 主程序 PROGRAM DEMO C WE COUNT FUNCTION AND JACOBIAN EVALUATIONS INTEGER NF NJ COMMON SCORE NF NJ C FOR PROBLEMS INTEGER N NB PARAMETER NB 100 DOUBLE PRECISION T H TOUT TREP Y NB B NB NB ATOLER RTOLER C WORKSPACE DECLARATION INTEGER NWORK PARAMETER NWORK 10000 INTEGER KOUNTI KOUNTD IWORK NWORK DOUBLE PRECISION DWORK NWORK C FOR BESIRK INTEGER INFO RESULT MAXIT ITER I STANDU DOUBLE PRECISION TREPRT HNEXT HMIN HMAX NSTEP ABSTOL NB RELTOL NB LOGICAL NUMJAC ALGEQN STANDT EXTERNAL STANDU STANDT C PROBLEM ROUTINE INTEGER REPT EXTERNAL FUNC REPT INTEGER M NS DOUBLE PRECISION AS 19 19 BS 19 19 XS 19 COMMON AB M NS AS BS XS C STEPSIZE IS ZERO BESIRK FINDS IT H 0 0D0 C INITIALIZE PROBLEM CALLLW WRITE MAIN WRITE AS 1 J J 1 NS M CALL INIT N T H TOUT TREP Y B NB NUMJAC ATOLER RTOLER C INITIALIZE COUNTERS NF 0 NJ 0 C SET INFO LEVEL INFO 0 C SET NEXT REPORT TIME TREPRT TREP C SET MINIMUM AND MAXIMUM STEPSIZE HMIN 1 0D 10 HMAX TOUT C SET NUMERICAL JACOBIAN STEP NSTEP 1 0D 3 C MAXIMUM NUMBER OF ITERATIONS MAXIT 500 C SET TOLERANCES FOR EACH VARIABLE DO I 1 N ABSTOL I ATOLER RELTOL I RTOLER ENDDO ALGEQN FALSE C SET WORKSPACE COUNTERS KOUNTI NWORK KOUNTD NWORK C SOLVE CALL INIT3 CALL BESIRK N INFO Y B NB T TOUT TREPRT H HNEXT HMIN HMAX RESULT NSTEP FUNC STANDU STANDT REPT NUMJAC MAXIT ITER ABSTOL RELTOL ALGEQN DWORK IWORK KOUNTD KOUNTI C REPORT WRITE WRITE BESIRK RETURNS RESULT WRITE FUNCTION EVALUATIONS NF WRITE JACOBIAN EVALUATIONS NJ WRITE C TERMINATE STOP END 子程序 SUBROUTINELW IMPLICIT REAL 8 A H O Z C EXTERNAL FEX JEX DIMENSION AS 19 19 BS 19 19 Q 19 19 XS 19 WS 19 C DIMENSION AU 19 19 BU 19 19 XA 19 WA 19 DIMENSION DIF1 19 DIF2 19 DIF3 19 ROOT 19 V1 19 V2 19 DIMENSION Y 99 ATOL 99 RWORK 10920 IWORK 120 COMMON AB N NS AS BS XS C COMMON BC Y1 YM2 YN1 C COMMON PARA PSAI SITA PE SAI PAK C N FOR SYMMETRIC COLLOCATION USED FOR PARTICLE AND C M FOR ASYMMETRIC COLLOCATION USED FOR COLUMN N 7 IW 1 IS 2 CALL COLL AS BS Q XS WS 19 N IW IS NS N 1 WRITE SYMMETRIC SITUATION WRITE POLYNOMIAL ROOTS WRITE XS I I 1 NS WRITE WRITE A MATRIX DO 20 I 1 NS 20 WRITE AS I J J 1 NS WRITE WRITE B MATRIX DO 30 I 1 NS 30 WRITE BS I J J 1 NS WRITE WRITE W MATRIX WRITE WS J J 1 NS WRITE ORXO WRITE AS 1 I I 1 NS N RETURN END SUBROUTINE FUNC N INFO F Y DFDY T CALCJ ERROR DW IW KOUNTD KOUNTI INTEGER N INFO ERROR KOUNTD KOUNTI IW DOUBLE PRECISION F N Y N DFDY N N T DW LOGICAL CALCJ INTEGER NF NJ COMMON SCORE NF NJ INTEGER M NS DOUBLE PRECISION AS 19 19 BS 19 19 XS 19 COMMON AB M NS AS BS XS NF NF 1 C WRITE FUNC C WRITE AS 1 J J 1 NS MM C PAUSE DO 99 J 1 M T1 0 0D0 DO 991 I 1 NS T1 T1 BS J I Y I 991CONTINUE T2 0 D0 DO 992 L 1 NS T2 T2 BS J L Y NS L 992CONTINUE F J T1 0 2 1 Y NS J EXP 20 20 Y J F NS J T2 0 3 1 Y NS J EXP 20 20 Y J 99 CONTINUE T0 0 D0 TT0 0 D0 DO 990 I 1 NS T0 T0 AS NS I Y I TT0 TT0 AS NS I Y NS I 990CONTINUE F NS Y NS 0 92 T0 F 2 NS TT0 RETURN END SUBROUTINE INIT N T H TOUT TREP Y B ND NUMJAC ATOLER RTOLER INTEGER N ND DOUBLE PRECISION T H TOUT TREP Y ND B ND ND ATOLER RTOLER LOGICAL NUMJAC INTEGER I DOUBLE PRECISION ZERO ONE PARAMETER ZERO 0 0D0 ONE 1 0D0 INTEGER M NS DOUBLE PRECISION AS 19 19 BS 19 19 XS 19 COMMON AB M NS AS BS XS N 2 NS T ZERO TOUT 1 0D0 TREP 0 1D0 DO I 1 NS Y I 1 0D0 Y NS I 0 0D0 ENDDO CALL MAT0 B ND N DO I 1 M B I I ONE B NS I NS I ONE ENDDO NUMJAC TRUE ATOLER 1 0D 6 RTOLER 0 0D0 RETURN END INTEGER FUNCTION REPT N ITER INFO T TREPRT Y E INTEGER N ITER INFO DOUBLE PRECISION T Y N E N TREPRT INTEGER M NS DOUBLE PRECISION AS 19 19 BS 19 19 XS 19 COMMON AB M NS AS BS XS OPEN 2 FILE JJM S DAE OUT WRITE 3 Z T WRITE 3 R WRITE 3 9 4X D11 5 XS I I 1 NS WRITE 3 T WRITE 3 9 4X D11 5 Y I I 1 NS WRITE 3 C WRITE 3 9 4X D11 5 Y NS I I 1 NS REPT 0 RETURN END C INCLUDE BESIRK FOR C INCLUDE MAIN FOR 输出结果 Z 0 000000000000000E 000 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 C 00000D 00 00000D 00 00000D 00 00000D 00 00000D 00 00000D 00 00000D 00 00000D 00 Z 7 437471009918388E 003 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 10015D 01 10015D 01 10015D 01 10015D 01 10010D 01 99895D 00 99577D 00 99382D 00 C 22625D 02 22625D 02 22625D 02 22615D 02 22460D 02 21962D 02 21537D 02 21463D 02 Z 9 055772950745512E 001 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 13939D 01 13771D 01 13500D 01 13175D 01 12848D 01 12568D 01 12372D 01 12284D 01 C 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 Z 1 000000000000000 R 14993D 00 33864D 00 51555D 00 67294D 00 80460D 00 90541D 00 97146D 00 10000D 01 T 13290D 01 13142D 01 12907D 01 12624D 01 12341D 01 12099D 01 11930D 01 11855D 01 C 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 二 利用非对称正交配置格式 1 不对称常微分方程程序 主程序 IMPLICIT REAL 8 A H O Z EXTERNAL FEX JEX C DIMENSION AS 19 19 BS 19 19 Q 19 19 XS 19 WS 19 DIMENSION AU 19 19 BU 19 19 XA 19 WA 19 DIMENSION DIF1 19 DIF2 19 DIF3 19 ROOT 19 V1 19 V2 19 DIMENSION Y 99 ATOL 99 RWORK 10920 IWORK 120 DOUBLE PRECISION Y1 Y2 YN1 YN2 COMMON AB M AU BU XA COMMON BC Y1 Y2 YN1 YN2 C N FOR SYMMETRIC COLLOCATION USED FOR PARTICLE AND C M FOR ASYMMETRIC COLLOCATION USED FOR COLUMN M 7 IW 1 C IS 2 C CALL COLL AS BS Q XS WS 19 N IW IS C NS N 1 C WRITE SYMMETRIC SITUATION CALL JCOBI 19 M 1 1 0 D0 0 D0 DIF1 DIF2 DIF3 ROOT NA M 2 WRITE ASYMMETRIC SITUATION WRITE THE POLYNOMIAL ROOTS DO 100 K 1 NA 100 XA K ROOT K WRITE XA I I 1 NA DO 140 I 1 NA CALL DFOPR 19 M 1 1 I 1 DIF1 DIF2 DIF3 ROOT V1 CALL DFOPR 19 M 1 1 I 2 DIF1 DIF2 DIF3 ROOT V2 WRITE DO 110 K 1 NA 110 AU I K V1 K DO 120 K 1 NA 120 BU I K V2 K WRITE AU I K K 1 NA WRITE BU I K K 1 NA 140 CONTINUE CALL RADAU 19 M 1 1 3 0 D0 0 D0 ROOT DIF1 V1 WRITE THE POLYNOMIAL WEIGHTS DO 150 K 1 NA 150 WA K V1 K WRITE WA I I 1 NA C CALCULATING THE PARAMETERS OF THE PROBLEM WHICH WILL BE USED C FOR THE DIMENSIONLESS FORM OF AND DEFINING OF THE PROBLEM NEQ 2 M LRW 22 9 NEQ NEQ 2 LIW 20 NEQ C INITIAL CONDITIONS DO 201 I 1 M Y I 1 D0 Y M I 0 D0 201 CONTINUE Y1 1 D0 YN1 1 D0 Y2 0 D0 YN2 0 D0 T 0 D0 DT 5 D 2 ITOL 2 RTOL 1 D 6 DO 203 I 1 NEQ ATOL I 1 D 6 203 CONTINUE ITASK 1 ISTATE 1 IOPT 0 MF 22 DO 240 IOUT 1 20 TOUT DT DFLOAT IOUT CALL LSODE FEX NEQ Y T TOUT ITOL RTOL ATOL ITASK ISTATE 1 IOPT RWORK LRW IWORK LIW JEX MF OPEN 2 FILE LW OUT WRITE 2 Z F8 4 T WRITE 2 R WRITE 2 10 4X D11 5 XA I I 1 NA WRITE 2 T WRITE 2 10 4X D11 5 Y1 Y I I 1 M YN1 WRITE 2 C C DO 205 I 1 M 2 WRITE 2 10 4X D11 5 Y2 Y M I I 1 M YN2 C WRITE 2 C 205 CONTINUE C 220 FORMAT 7H AT T D12 4 6H Y 3D15 7 IF ISTATE LT 0 GO TO 280 240 CONTINUE WRITE 3 260 IWORK 11 IWORK 12 IWORK 13 260 FORMAT 12H NO STEPS I4 11H NO F S I4 11H NO J S I4 STOP 280 WRITE 3 290 ISTATE 290 FORMAT 22H ERROR HALT ISTATE I3 STOP END 子程序 SUBROUTINE FEX NEQ T Y YP C VARIABLES IN THE ORDER OF Y 1 TO Y M FOR THE COLUMN COLLOCATION C POINTS OF 2 TO M 1 Y M 1 TO Y M N FOR PARTICLE COLLOCATION POINTS C FROM 1 TO N IN COLUMN COLLOCATION POINT 1 AND Y M N M 1 1 TO C Y M N M 2 FOR PARTICLE COLLOCATION POINTS FROM 1 TO N IN COLUMN C COLLOCATION POINT M 2 IMPLICIT REAL 8 A H O Z DIMENSION Y 19 19 YP 19 19 AU 19 19 BU 19 19 XA 19 DOUBLE PRECISION Y1 Y2 YN1 YN2 COMMON AB M AU BU XA COMMON BC Y1 Y2 YN1 YN2 T0 0 D0 TT0 0 D0 TN1 0 D0 TTN1 0 D0 DO 990 I 2 M 1 T0 T0 AU 1 I Y I 1 TN1 TN1 AU M 2 I Y I 1 TT0 TT0 AU 1 I Y M I 1 TTN1 TTN1 AU M 2 I Y M I 1 990CONTINUE Y1 0 92 TN1 1 AU M 2 M 2 T0 AU 1 M 2 AU M 2 1 AU M 2 M 2 1 AU 1 1 AU 1 M 2 YN1 T0 AU 1 1 Y1 AU 1 M 2 Y2 AU 1 M 2 TTN1 AU M 2 M 2 TT0 AU 1 1 AU 1 M 2 AU M 2 1 AU M 2 M 2 YN2 TTN1 AU M 2 1 Y2 AU M 2 M 2 C DO 99 J 2 M 1 T1 AU J 1 Y1 AU J M 2 YN1 T2 AU J 1 Y2 AU J M 2 YN2 TB1 BU J 1 Y1 BU J M 2 YN1 TB2 BU J 1 Y2 BU J M 2 YN2 DO 991 I 1 M T1 T1 AU J I 1 Y I T2 T2 AU J I 1 Y M I TB1 TB1 BU J I 1 Y I TB2 TB2 BU J I 1 Y M I 991CONTINUE YP J 1 1 XA J T1 TB1 0 2 1 Y M J 1 EXP 20 20 Y J 1 YP M J 1 1 XA J T2 TB2 0 3 1 Y M J 1 EXP 20 20 Y J 1 99 CONTINUE RETURN END SUBROUTINE JEX NEQ T Y ML MU PD NRPD DOUBLE PRECISION PD T Y DIMENSION Y NEQ PD NRPD NEQ RETURN END 输出结果 Z 0500 R 00000D 00 25446D 01 12923D 00 29708D 00 50000D 00 70292D 00 87077D 00 97455D 00 10000D 01 T 10113D 01 10109D 01 10109D 01 10105D 01 10090D 01 10047D 01 99758D 00 99133D 00 98957D 00 C 17090D 01 16543D 01 16530D 01 16448D 01 16127D 01 15385D 01 14587D 01 14286D 01 14705D 01 Z 1000 R 00000D 00 25446D 01 12923D 00 29708D 00 50000D 00 70292D 00 87077D 00 97455D 00 10000D 01 T 10221D 01 10220D 01 10216D 01 10199D 01 10154D 01 10077D 01 99859D 00 99184D 00 99010D 00 C 36197D 01 35926D 01 35774D 01 35082D 01 33496D 01 31339D 01 29744D 01 29256D 01 29436D 01 Z 9500 R 00000D 00 25446D 01 12923D 00 29708D 00 50000D 00 70292D 00 87077D 00 97455D 00 10000D 01 T 13564D 01 13661D 01 13633D 01 13510D 01 13236D 01 12839D 01 12430D 01 12147D 01 12010D 01 C 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 Z 1 0000 R 00000D 00 25446D 01 12923D 00 29708D 00 50000D 00 70292D 00 87077D 00 97455D 00 10000D 01 T 13234D 01 13325D 01 13299D 01 13184D 01 12931D 01 12563D 01 12184D 01 11922D 01 11796D 01 C 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 10000D 01 NO STEPS 227 NO F S 700 NO J S 29 2 不对称代数微分方程程序 主程序为 PROGRAM DEMO C WE COUNT FUNCTION AND JACOBIAN EVALUATIONS INTEGER NF NJ COMMON SCORE NF NJ C FOR PROBLEMS INTEGER N NB PARAMETER NB 100 DOUBLE PRECISION T H TOUT TREP Y NB B NB NB ATOLER RTOLER C WORKSPACE DECLARATION INTEGER NWORK PARAMETER NWORK 10000 INTEGER KOUNTI KOUNTD IWORK NWORK DOUBLE PRECISION DWORK NWORK C FOR BESIRK INTEGER INFO RESULT MAXIT ITER I STANDU DOUBLE PRECISION TREPRT HNEXT HMIN HMAX NSTEP ABSTOL NB RELTOL NB LOGICAL NUMJAC ALGEQN STANDT EXTERNAL STANDU STANDT C PROBLEM ROUTINE INTEGER REPT EXTERNAL FUNC REPT INTEGER M NA DOUBLE PRECISION AU 19 19 BU 19 19 XA 19 COMMON AB M NA AU BU XA C STEPSIZE IS ZERO BESIRK FINDS IT H 0 0D0 C INITIALIZE PROBLEM CALL LW WRITE MAIN WRITE AU 1 I I 1 NA M NA CALL INIT N T H TOUT TREP Y B NB NUMJAC ATOLER RTOLER C INITIALIZE COUNTERS NF 0 NJ 0 C SET INFO LEVEL INFO 0 C SET NEXT REPORT TIME TREPRT TREP C SET MINIMUM AND MAXIMUM STEPSIZE HMIN 1 0D 10 HMAX TOUT C SET NUMERICAL JACOBIAN STEP NSTEP 1 0D 3 C MAXIMUM NUMBER OF ITERATIONS MAXIT 500 C SET TOLERANCES FOR EACH VARIABLE DO I 1 N ABSTOL I ATOLER RELTOL I RTOLER ENDDO ALGEQN FALSE C SET WORKSPACE COUNTERS KOUNTI NWORK KOUNTD NWORK C SOLVE CALL INIT3 CALL BESIRK N INFO Y B NB T TOUT TREPRT H HNEXT HMIN HMAX RESULT NSTEP FUNC STANDU STANDT REPT NUMJAC MAXIT ITER ABSTOL RELTOL ALGEQN DWORK IWORK KOUNTD KOUNTI C REPORT WRITE WRITE BESIRK RETURNS RESULT WRITE FUNCTION EVALUATIONS NF WRITE JACOBIAN EVALUATIONS NJ WRITE C TERMINATE STOP END 子程序为 SUBROUTINE FUNC N INFO F Y DFDY T CALCJ ERROR DW IW KOUNTD KOUNTI INTEGER M NA DOUBLE PRECISION AU 19 19 BU 19 19 XA 19 COMMON AB M NA AU BU XA INTEGER N INFO ERROR KOUNTD KOUNTI IW DOUBLE PRECISION F N Y N DFDY N N T DW LOGICAL CALCJ INTEGER NF NJ COMMON SCORE NF NJ NF NF 1 C WRITE AU 1 I I 1 NA M NA DO 99 J 2 M 1 T1 0 0D0 T2 0 D0 T3 0 D0 T4 0 D0 DO 991 I 1 NA T1 T1 AU J I Y I T2 T2 BU J I Y I T3 T3 AU J I Y NA I T4 T4 BU J I Y NA I 991 CONTINUE F J 1 XA J T1 T2 0 2 1 Y NA J EXP 20 20 Y J F NA J 1 XA J T3 T4 0 3 1 Y NA J EXP 20 20 Y J 99CONTINUE T0 0 D0 TT0 0 D0 TNAA 0 D0 TTNAA 0 D0 DO 992 I 1 NA T0 T0 AU 1 I Y I TT0 TT0 AU 1 I Y NA I TNAA TNAA AU NA I Y I TTNAA TTNAA AU NA I Y NA I 992 CONTINUE F 1 T0 F NA TNAA Y NA 0 92 F NA 1 TT0 F 2 NA TTNAA RETURN END SUBROUTINE INIT N T H TOUT TREP Y B ND NUMJAC ATOLER RTOLER INTEGER N ND DOUBLE PRECISION T H TOUT TREP Y ND B ND ND ATOLER RTOLER LOGICAL NUMJAC INTEGER I DOUBLE PRECISION ZERO ONE PARAMETER ZERO 0 0D0 ONE 1 0D0 INTEGER M NA DOUBLE PRECISION AU 19 19 BU 19 19 XA 19 COMMON AB M NA AU BU XA N 2 NA T ZERO TOUT 1 0D0 TREP 0 1D0 DO I 1 NA Y I 1 D0 Y NA I 0 D0 ENDDO CALL MAT0 B ND N DO I 2 M 1 B I I ONE B NA I NA I ONE ENDDO NUMJAC TRUE ATOLER 1 0D 6 RTOLER 0 0D0 RETURN END INTEGER FUNCTION REPT N ITER INFO T TREPRT Y E INTEGER N ITER INFO DOUBLE PRECISION T Y N E N TREPRT INTEGER M NA DOUBLE PRECISION AU 19 19 BU 19 19 XA 19 COMMON AB M NA AU BU XA OPEN FILE JJM AS DAE OUT WRITE 3 Z T WRITE 3 R WRITE 3 10 4X D11 5 XA I I 1 NA WRITE 3 T WRITE 3 10 4X D11 5 Y I I 1 NA WRITE 3 C WRITE 3 10 4X D11 5 Y NA I I 1 NA REPT 0 RETURN END SUBROUTINELW IMPLICIT REAL 8 A H O Z EXTERNAL FEX JEX C DIMENSION AS 19 19 BS 19 19 Q 19 19 XS 19 WS 19 DIMENSION AU 19 19 BU 19 19 XA 19 WA 19 DIMENSION DIF1 19 DIF2 19 DIF3 19 ROOT 19 V1 19 V2 19 DIMENSION Y 99 ATOL 99 RWORK 10920 IWORK 120 COMMON AB M NA AU BU XA C COMMON BC Y1 YM2 YN1 C COMMON PARA PSAI SITA PE SAI PAK C N FOR SYMMETRIC COLLOCATION USED FOR PARTICLE AND C M