(original) (raw)

SUBROUTINE RADAUP(N,FCN,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JAC ,IJAC,MLJAC,MUJAC, & MAS ,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C ---------------------------------------------------------- C NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) C SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS C M*Y'=F(X,Y). C THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I) C OR EXPLICIT (M=I). C THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA) C OF ORDER 5, 9 0R 13 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT. C C.F. SECTION IV.8 C C AUTHORS: E. HAIRER AND G. WANNER C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES C CH-1211 GENEVE 24, SWITZERLAND C E-MAIL: Ernst.Hairer@math.unige.ch C Gerhard.Wanner@math.unige.ch C C THIS CODE IS PART OF THE BOOK: C E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL C EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, C SPRINGER-VERLAG 1991, SECOND EDITION 1996. C C VERSION OF JULY 9, 1996 C C INPUT PARAMETERS C ---------------- C N DIMENSION OF THE SYSTEM C C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE C VALUE OF F(X,Y): C SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR) C DOUBLE PRECISION X,Y(N),F(N) C F(1)=... ETC. C RPAR, IPAR (SEE BELOW) C C X INITIAL X-VALUE C C Y(N) INITIAL VALUES FOR Y C C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) C C H INITIAL STEP SIZE GUESS; C FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, C H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD. C THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS C QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6). C C RTOL,ATOL RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY C CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. C C ITOL SWITCH FOR RTOL AND ATOL: C ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. C THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF C Y(I) BELOW RTOL*ABS(Y(I))+ATOL C ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. C THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW C RTOL(I)*ABS(Y(I))+ATOL(I). C C JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES C THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y C (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY C A DUMMY SUBROUTINE IN THE CASE IJAC=0). C FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM C SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR) C DOUBLE PRECISION X,Y(N),DFY(LDFY,N) C DFY(1,1)= ... C LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS C FURNISHED BY THE CALLING PROGRAM. C IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO C BE FULL AND THE PARTIAL DERIVATIVES ARE C STORED IN DFY AS C DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) C ELSE, THE JACOBIAN IS TAKEN AS BANDED AND C THE PARTIAL DERIVATIVES ARE STORED C DIAGONAL-WISE AS C DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). C C IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN: C IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE C DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. C IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC. C C MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: C MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. C 0<=MLJAC= NUMBER OF NON-ZERO DIAGONALS BELOW C THE MAIN DIAGONAL). C C MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON- C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). C NEED NOT BE DEFINED IF MLJAC=N. C C ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- C ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - C C MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- C MATRIX M. C IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY C MATRIX AND NEEDS NOT TO BE DEFINED; C SUPPLY A DUMMY SUBROUTINE IN THIS CASE. C IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM C SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR) C DOUBLE PRECISION AM(LMAS,N) C AM(1,1)= .... C IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED C AS FULL MATRIX LIKE C AM(I,J) = M(I,J) C ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED C DIAGONAL-WISE AS C AM(I-J+MUMAS+1,J) = M(I,J). C C IMAS GIVES INFORMATION ON THE MASS-MATRIX: C IMAS=0: M IS SUPPOSED TO BE THE IDENTITY C MATRIX, MAS IS NEVER CALLED. C IMAS=1: MASS-MATRIX IS SUPPLIED. C C MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: C MLMAS=N: THE FULL MATRIX CASE. THE LINEAR C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. C 0<=MLMAS= NUMBER OF NON-ZERO DIAGONALS BELOW C THE MAIN DIAGONAL). C MLMAS IS SUPPOSED TO BE .LE. MLJAC. C C MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). C NEED NOT BE DEFINED IF MLMAS=N. C MUMAS IS SUPPOSED TO BE .LE. MUJAC. C C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE C NUMERICAL SOLUTION DURING INTEGRATION. C IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. C IT MUST HAVE THE FORM C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N, C RPAR,IPAR,IRTRN) C DOUBLE PRECISION X,Y(N),CONT(LRC) C .... C SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS C THE FIRST GRID-POINT). C "XOLD" IS THE PRECEEDING GRID-POINT. C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN C IS SET <0, RADAU5 RETURNS TO THE CALLING PROGRAM. C C ----- CONTINUOUS OUTPUT: ----- C DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION C FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH C THE FUNCTION C >>> CONTRP(I,S,CONT,LRC) <<< C WHICH PROVIDES AN APPROXIMATION TO THE I-TH C COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE C S SHOULD LIE IN THE INTERVAL [XOLD,X]. C C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLOUT: C IOUT=0: SUBROUTINE IS NEVER CALLED C IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT. C C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". C WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS C FOR THE CODE. FOR STANDARD USE OF THE CODE C WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE C CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE. C WORK(8),..,WORK(LWORK) SERVE AS WORKING SPACE C FOR ALL VECTORS AND MATRICES. C "LWORK" MUST BE AT LEAST C N*(LJAC+LMAS+NS*LE+3*NS+3)+20 C WHERE C NS=IWORK(11) (SEE BELOW) C AND C LJAC=N IF MLJAC=N (FULL JACOBIAN) C LJAC=MLJAC+MUJAC+1 IF MLJAC0 THEN "LWORK" MUST BE AT LEAST C N*(LJAC+3*NS+3)+(N-M1)*(LMAS+NS*LE)+20 C WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE C NUMBER N CAN BE REPLACED BY N-M1. C C LWORK DECLARED LENGTH OF ARRAY "WORK". C C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". C IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS C FOR THE CODE. FOR STANDARD USE, SET IWORK(1),.., C IWORK(20) TO ZERO BEFORE CALLING. C IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA. C "LIWORK" MUST BE AT LEAST C (2+(NS-1)/2)*N+20. C C LIWORK DECLARED LENGTH OF ARRAY "IWORK". C C RPAR, IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH C CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING C PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. C C ---------------------------------------------------------------------- C C SOPHISTICATED SETTING OF PARAMETERS C ----------------------------------- C SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK C WELL. THEY MAY BE DEFINED BY SETTING WORK(1),... C AS WELL AS IWORK(1),... DIFFERENT FROM ZERO. C FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: C C IWORK(1) IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN C MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY C ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN. C IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<n) 3="" c="" and="" not="" for="" implicit="" systems="" (imas="1)." iwork(2)="" this="" is="" the="" maximal="" number="" of="" allowed="" steps.="" default="" value="" (for="" 100000.="" iwork(3)="" maximum="" newton="" iterations="" solution="" system="" in="" each="" step.="" 7+(ns-3)*2.="" ns="" stages="" (see="" iwork(11)).="" iwork(4)="" if="" iwork(4).eq.0="" extrapolated="" collocation="" taken="" as="" starting="" newton's="" method.="" iwork(4).ne.0="" zero="" values="" are="" used.="" latter="" recommended="" method="" has="" difficulties="" with="" convergence="" (this="" case="" when="" nstep="" larger="" than="" naccpt="" +="" nrejct;="" see="" output="" param.).="" following="" parameters="" important="" differential-algebraic="" index=""> 1. C THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT C THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER. C IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE C MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2. C C IWORK(5) DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR C ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM. C DEFAULT IWORK(5)=N. C C IWORK(6) DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0. C C IWORK(7) DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0. C C IWORK(8) SWITCH FOR STEP SIZE STRATEGY C IF IWORK(8).EQ.1 MOD. PREDICTIVE CONTROLLER (GUSTAFSSON) C IF IWORK(8).EQ.2 CLASSICAL STEP SIZE CONTROL C THE DEFAULT VALUE (FOR IWORK(8)=0) IS IWORK(8)=1. C THE CHOICE IWORK(8).EQ.1 SEEMS TO PRODUCE SAFER RESULTS; C FOR SIMPLE PROBLEMS, THE CHOICE IWORK(8).EQ.2 PRODUCES C OFTEN SLIGHTLY FASTER RUNS C C IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT C Y(I)' = Y(I+M2) FOR I=1,...,M1, C WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME C CAN BE ACHIEVED BY SETTING THE PARAMETERS IWORK(9) AND IWORK(10). C E.G., FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE C VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2. C FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS: C - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE C JACOBIAN HAVE TO BE STORED C IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL C DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J) C FOR I=1,N-M1 AND J=1,N. C ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM ) C DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2) C FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM. C - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL C 0<=MLJAC</n)>