SUBROUTINE INITAL C... C... KROGH'S ODE PROBLEM C... C... KROGH'S PROBLEM IS A NONLINEAR ODE PROBLEM WITH AN EXACT C... SOLUTION. THE NUMBER OF ODES AND THE STIFFNESS OF THESE C... EQUATIONS CAN BE VARIED BY THE CHOICE OF PARAMETERS. IN C... THIS CASE, N = 20, AND THE STIFFNESS IS SET BY CONSTANTS C... READ FROM A DATA FILE. C... C... THE NONLINEAR ODE PROBLEM PROPOSED BY KROGH (GEAR (1), PP. C... 218-220), IS AN IDEAL TEST PROBLEM FOR ODE SOFTWARE FOR THE C... FOLLOWING REASONS: C... C... (1) THE ORDER OF THE PROBLEM (NUMBER OF ODES) CAN BE C... SELECTED BY THE USER (IN THE PRESENT CASE, N = 20). C... C... (2) THE STIFFNESS OF THE ODES CAN BE SELECTED BY THE USER. C... THIS POINT REQUIRES SOME CLARIFICATION SINCE THE ODES C... ARE NONLINEAR, BUT STIFFNESS IS A LINEAR CONCEPT BASED C... ON EIGENVALUES (A STIFF SYSTEM IS GENERALLY DEFINED AS: C... C... (PROBLEM TIME)*(LARGEST EIGENVALUE) C... C... IS MUCH GREATER THAN 1) C... C... GENERALLY, THE CONCEPT OF STIFFNESS APPLIED TO NONLINEAR C... PROBLEMS MEANS APPLIED TO THE LINEARIZED VERSION OF THE C... NONLINEAR EQUATIONS, FOR WHICH EIGENVALUES CAN BE COM- C... PUTED. HOWEVER, SINCE LINEARIZATION IS PERFORMED AROUND C... A BASE POINT, THE EIGENVALUES OF THE LINEARIZED SYSTEM C... WILL CHANGE WITH THE BASE POINT, I.E., THE LOCAL EIGEN- C... VALUES OF A NONLINEAR SYSTEM ARE NOT CONSTANT. C... C... THE NONLINEAR ODES IN THE KROGH PROBLEM HAVE IN THEIR C... ANALYTICAL SOLUTION EXPONENTIALS OF THE FORM EXP(B(I)*T) C... WHERE T IS THE INDEPENDENT VARIABLE AND THE B(I), I = C... 1, 2,..., N ARE PARAMETERS SELECTED BY THE USER. THUS C... THE B(I) ARE ANALOGOUS TO EIGENVALUES, BUT THE EXPONENT- C... IALS (EXP(B(I)*T)) APPEAR NONLINEARLY IN THE ANALYTICAL C... SOLUTION (RATHER THAN AS A LINEAR COMBINATION IN THE C... ANALYTICAL SOLUTION OF A SYSTEM OF CONSTANT COEFFICIENT C... LINEAR ODES). C... C... IN THIS SENSE, THE STIFFNESS OF THE KROGH TEST PROBLEM C... CAN BE SPECIFIED BY THE USER BY SELECTING THE B(I). C... C... (3) AS THE PRECEDING DISCUSSION INDICATES, THE KROGH TEST C... HAS AN ANALYTICAL (EXACT) SOLUTION, WHICH CAN BE USED C... AS A STANDARD IN EVALUATING A NUMERICAL INTEGRATOR, I.E., C... THE ERRORS IN NUMERICAL SOLUTIONS CAN BE COMPUTED. C... C... ALL OF THE FEATURES OF THE KROGH TEST PROBLEM ARE ILLUSTRATED IN C... THE FOLLOWING CODING FOR TWO NONSTIFF CASES: C... C... (1) WE TAKE B(I) = -I, I = 1, 2,..., 20 C... C... THESE CONSTANTS ARE READ FROM A DATA FILE IN SUBROUTINE C... INITAL AT THE BEGINNING OF THE NUMERICAL SOLUTION, AND C... THE INTEGRATION IS PERFORMED BY BY THE SECOND-ORDER C... ADAMS-BASHFORTH METHOD C... C... (1) WE TAKE B(I) = I, I = 1, 2,..., 20 (NOTE THE ONLY CHANGE C... IS IN THE SIGN OF B(I)). C... C... THE MATHEMATICAL DETAILS OF THE PROBLEM ARE GIVEN IN REFERENCE C... (1) CITED PREVIOUSLY. THEY CAN ALSO BE INFERRED FROM THE C... FOLLOWING CODING. C... COMMON /T/ T COMMON /Y/ Y(20) COMMON /F/ DYDT(20) COMMON /C/ BETA(20), U(20,20), B(20,20), N COMMON /IO/ NI, NO C... C... READ THE NUMBER OF ODES AND THE CONSTANT ARRAY, BETA, WHICH C... DEFINES THE PROBLEM TIME SCALE AND STIFFNESS READ(NI,100)N DO 10 J=1,N READ(NI,101)BETA(J) 10 CONTINUE WRITE(NO,102)N,(BETA(J),J=1,N) C... C... INITIAL CONDITIONS DO 20 J=1,N Y(J) = -1.0 20 CONTINUE C... C... SET THE CONSTANT ARRAY, U TERM = 2.0/FLOAT(N) DO 30 I=1,N DO 40 J=1,N U(I,J) = TERM 40 CONTINUE U(I,I) = TERM - 1.0 30 CONTINUE C... C... SET THE CONSTANT ARRAY, B DO 50 J=1,N DO 60 I=1,N SUM = 0.0 DO 70 K=1,N SUM = SUM + U(I,K)*BETA(K)*U(K,J) 70 CONTINUE B(I,J) = SUM 60 CONTINUE 50 CONTINUE C... C... FORMATS ARE LISTED HERE 100 FORMAT(I5) 101 FORMAT(F20.5) 102 FORMAT(/,11X,'N = ',I3,10X,'BETA',20(/,17X,F15.3)) RETURN END SUBROUTINE DERV COMMON /T/ T COMMON /Y/ Y(20) COMMON /F/ DYDT(20) COMMON /C/ BETA(20), U(20,20), B(20,20), N DIMENSION W(20) C... C... CALCULATE THE ARRAY, W DO 10 I=1,N SUM = 0.0 DO 20 J=1,N SUM = SUM + U(I,J)*Y(J) 20 CONTINUE W(I) = SUM 10 CONTINUE C... C... CALCULATE THE N ODE DERIVATIVES DO 30 I=1,N SUM = 0.0 DO 40 J=1,N SUM = SUM + U(I,J)*W(J)*W(J) - B(I,J)*Y(J) 40 CONTINUE DYDT(I) = SUM 30 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) COMMON /T/ T COMMON /Y/ Y(20) COMMON /F/ DYDT(20) COMMON /C/ BETA(20), U(20,20), B(20,20), N DIMENSION W(20), YEX(20) C... C... ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY SUBROUTINE JMAP C... (A, X, XOLD, F, FOLD) AND SUBROUTINE EIGEN (WI, WR, Z, RW, IW) REAL 1 A(20,20), WR(20), WI(20), Z(20,20), RW(20), 2 X(20), XOLD(20), 3 F(20), FOLD(20) INTEGER IW(20) C... C... THE FOLLOWING EQUIVALENCE FACILITATES TRANSFER OF THE DEPENDENT C... VARIABLE VECTOR IN /Y/ TO ARRAY Y, AND THE DERIVATIVE VECTOR IN C... IN /F/ TO ARRAY F FOR USE IN THE SUBSEQUENT CALL TO SUBROUTINE C... JMAP EQUIVALENCE (Y(1),X(1)),(DYDT(1),F(1)) C... C... INITIALIZE A COUNTER FOR THE CALLS TO SUBROUTINES JMAP AND EIGEN DATA IP/0/ IP=IP+1 C... C... MONITOR THE PROGRESS OF THE SOLUTION ON THE SCREEN WRITE(*,*)' IP = ',IP C... C... MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE C... DERV, AND COMPUTE ITS TEMPORAL EIGENVALUES. SUBROUTINE EIGEN C... CALLS A SERIES OF EISPACK ROUTINES TO COMPUTE THE TEMPORAL C... EIGENVALUES, AND OPTIONALLY, THE EIGENVECTORS OF THE 20-ODE C... SYSTEM JACOBIAN MATRIX (THE CALLS TO JMAP AND EIGEN ARE MADE C... ONLY AT THE BEGINNING AND END OF EACH RUN TO REDUCE THE VOLUME C... OF OUTPUT) IF(IP.EQ.1)THEN CALL JMAP(N,A, X,XOLD,F,FOLD ) CALL EIGEN(N,A,WR, WI,Z, RW,IW) C... C... PRINT A HEADING FOR THE SOLUTION WRITE(NO,2) 2 FORMAT(//,3X,'TIME',5X,'Y(1)',3X,'YEX(1)',5X,'PER1', 1 4X,'Y(20)',2X,'YEX(20)',4X,'PER20', 2 3X,'PERMAX') END IF C... C... CALCULATE THE N EXACT SOLUTIONS DO 10 J=1,N BETAT = BETA(J)*T IF (BETAT.GT.200.) THEN W(J) = 0.0 ELSE IF (BETAT.LT.-200.) THEN W(J) = BETA(J) ELSE W(J) = BETA(J)/(1.0 - (1.0 + BETA(J))*EXP(BETAT)) END IF 10 CONTINUE DO 20 I=1,N SUM = 0.0 DO 30 J=1,N SUM = SUM + U(I,J)*W(J) 30 CONTINUE YEX(I) = SUM 20 CONTINUE C... C... CALCULATE THE MAXIMUM ERROR ERRMAX = ABS((Y(1) - YEX(1))/YEX(1))*100. DO 40 J=2,N ERR = ABS((Y(J) - YEX(J))/YEX(J))*100. IF (ERR.GT.ERRMAX) ERRMAX = ERR 40 CONTINUE C... C... PRINT THE NUMERICAL AND EXACT SOLUTIONS FOR I = 1 AND N, AND C... THE CORRESPONDING ERRORS, AND THE MAXIMUM ERROR PER1=ABS((Y(1)-YEX(1))/YEX(1))*100. PERN=ABS((Y(N)-YEX(N))/YEX(N))*100. PERMAX=ERRMAX WRITE(NO,3)T,Y(1),YEX(1),PER1,Y(N),YEX(N),PERN,PERMAX 3 FORMAT(F7.2,7F9.4) C... C... MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE C... DERV, AND COMPUTE ITS TEMPORAL EIGENVALUES AT THE END OF THE RUN IF(IP.EQ.11)THEN CALL JMAP(N,A, X,XOLD,F,FOLD ) CALL EIGEN(N,A,WR, WI,Z, RW,IW) END IF C... C... RESET IP AT THE END OF EACH RUN IF(IP.EQ.11)IP=0 RETURN END