!**************************************************************** !* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY * !* DIFFERENTIAL EQUATIONS Y'=F(X,Y) BY ROSENBROCK METHOD. * !* ------------------------------------------------------------ * !* SAMPLE RUN: * !* Example #1: * !* (Solve set of differential equations (N=2): * !* F(1) = Y(1) * Y(2) + COS(X) - HALF * SIN(TWO * X) * !* F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + SIN(X)) * !* Find values of F(1), F(2) at X=1.5). * !* * !* SOLUTION AT X= 1.50000000000000 * !* Y(1) = 0.12359935E+01 * !* Y(2) = -0.10494372E+00 * !* * !* LAST STEP SIZE = 4.150113101356574E-002 * !* ERROR CODE = 1 * !* * !* Example #2: * !* (Solve set of differential equations (N=5): * !* F(1) = Y(2) * !* F(2) = Y(3) * !* F(3) = Y(4) * !* F(4) = Y(5) * !* F(5) = (45.d0 * Y(3) * Y(4) * Y(5) - * !* 40.d0 * Y(4) * Y(4) * Y(4)) / (NINE * Y(3) * Y(3)) * !* Find values of F(1), F(2), ..., F(5) at X=1.5). * !* * !* SOLUTION AT X= 1.50000000000000 * !* Y(1) = 0.43639610E+01 * !* Y(2) = 0.40000000E+01 * !* Y(3) = 0.28284271E+01 * !* Y(4) = 0.14790900E-10 * !* Y(5) = -0.37712362E+01 * !* * !* LAST STEP SIZE = 3.825256949194526E-003 * !* ERROR CODE = 1 * !* ------------------------------------------------------------ * !* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * !* [BIBLI 18]. * !* * !* F90 Release 1.0 By J-P Moreau, Paris * !* (www.jpmoreau.fr) * !**************************************************************** ! LIST OF USED SUBROUTINES (HERE INCLUDED): ! ======================================== ! ROS4, RO4COR, SHAMP, GRK4A, GRK4T, VELDD, VELDS, LSTAB, ! DEC, DECB, SOL, SOLB. ! (Other used subroutines/functions are standard Fortran Library). ! --------------------------------------------------------------- PROGRAM TROS4 REAL*8 X,XEND, H,RTOL,ATOL REAL*8,POINTER :: Y(:) REAL*8,POINTER :: WORK(:) INTEGER,POINTER :: IWORK(:) !Initialize parameters (see ROS4) N=2 !DIMENSION OF THE SYSTEM (N=5 for #2) IFCN=1 !FCN(N,X,Y,F) MAY DEPEND ON X X=0.d0 !INITIAL X-VALUE XEND=1.5d0 !FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) H=0.001d0 !INITIAL STEP SIZE GUESS RTOL=1.d-6 !RELATIVE ERROR TOLERANCE (HERE SCALAR) ATOL=1.d-8 !ABSOLUTE ERROR TOLERANCE (HERE SCALAR) !1.d-10 for both in #2. ITOL=0 !BOTH RTOL AND ATOL ARE SCALARS IJAC=0 !JACOBIAN IS COMPUTED INTERNALLY BY FINITE !DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED MLJAC=N !JACOBIAN IS A FULL MATRIX. THE LINEAR ALGEBRA !IS DONE BY FULL-MATRIX GAUSS-ELIMINATION IDFX=0 !DF/DX IS COMPUTED INTERNALLY BY FINITE !DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED IMAS=0 !M IS SUPPOSED TO BE THE IDENTITY !MATRIX, MAS IS NEVER CALLED MLMAS=N !MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ALGEBRA !IS DONE BY FULL-MATRIX GAUSS-ELIMINATION IOUT=0 !SUBROUTINE SOLOUT IS NEVER CALLED LE1=N !IF MLJAC=N (FULL JACOBIAN) LJAC=N !IF MLJAC=N (FULL JACOBIAN) LMAS=0 !IF IMAS=0 LIWORK= N+2 !DECLARED LENGTH OF ARRAY "IWORK" LWORK = N*(LJAC+LMAS+LE1+8)+5 !DECLARED LENGTH OF ARRAY "LWORK" !dynamic allocations allocate(Y(1:N),stat=ialloc) allocate(WORK(1:LWORK),stat=ialloc) allocate(IWORK(1:LIWORK),stat=ialloc) WORK=0.d0 !This triggers default values (see ROS4) IWORK=0 Y(1)=0.5d0 !INITIAL VALUES FOR Y Y(2)=0.5d0 !In #2, Y(1) = Y(2) = ... = Y(5) = 1.d0 !call Rosenbrock subroutine with appropriate parameters !(here, FCN has been removed from input parameters). CALL ROS4(N,IFCN,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX, & MAS ,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,IDID) !print results print *,' ' print *,' SOLUTION AT X=', X do I=1,N write(*,10) I, Y(I) end do print *,' ' print *,' LAST STEP SIZE =', H print *,' ERROR CODE =', IDID print *,' ' 10 format(' Y(',I1,') = ',E15.8) END !of main program !define example #1 SUBROUTINE FCN(N,X,Y,F) parameter(HALF=0.5d0,ONE=1.d0,TWO=2.d0) REAL*8 X,Y(N),F(N) F(1) = Y(1) * Y(2) + DCOS(X) - HALF * DSIN(TWO * X) F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + DSIN(X)) RETURN END !define example #2 !SUBROUTINE FCN(N,X,Y,F) !parameter(NINE=9.d0) !REAL*8 X,Y(N),F(N) ! F(1) = Y(2) ! F(2) = Y(3) ! F(3) = Y(4) ! F(4) = Y(5) ! F(5) = (45.d0 * Y(3) * Y(4) * Y(5) - & ! 40.d0 * Y(4) * Y(4) * Y(4)) / (NINE * Y(3) * Y(3)) !RETURN !END !********************************************************************** SUBROUTINE ROS4(N,IFCN,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX, & MAS ,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,IDID) ! --------------------------------------------------------------------- ! NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) ! SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(X,Y). ! THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4 ! (WITH STEP SIZE CONTROL). ! C.F. SECTION IV.7 ! ! AUTHORS: E. HAIRER AND G. WANNER ! UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES ! CH-1211 GENEVE 24, SWITZERLAND ! E-MAIL: HAIRER@CGEUGE51.BITNET, WANNER@CGEUGE51.BITNET ! ! THIS CODE IS PART OF THE BOOK: ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, ! SPRINGER-VERLAG (1990) ! ! VERSION OF OCTOBER 12, 1990 ! ! INPUT PARAMETERS ! ---------------- ! N DIMENSION OF THE SYSTEM ! ! FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE ! VALUE OF F(X,Y): ! SUBROUTINE FCN(N,X,Y,F) ! REAL*8 X,Y(N),F(N) ! F(1)=... ETC. ! ! IFCN GIVES INFORMATION ON FCN: ! IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS) ! IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS) ! ! X INITIAL X-VALUE ! ! Y(N) INITIAL VALUES FOR Y ! ! XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) ! ! H INITIAL STEP SIZE GUESS; ! FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, ! H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD. ! THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY ! ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW ! STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE. ! (IF H=0.D0, THE CODE PUTS H=1.D-6). ! ! RTOL,ATOL RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY ! CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. ! ! ITOL SWITCH FOR RTOL AND ATOL: ! ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. ! THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF ! Y(I) BELOW RTOL*ABS(Y(I))+ATOL ! ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. ! THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW ! RTOL(I)*ABS(Y(I))+ATOL(I). ! ! JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES ! THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y ! (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY ! A DUMMY SUBROUTINE IN THE CASE IJAC=0). ! FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM: ! SUBROUTINE JAC(N,X,Y,DFY,LDFY) ! REAL*8 X,Y(N),DFY(LDFY,N) ! DFY(1,1)= ... ! LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS ! FURNISHED BY THE CALLING PROGRAM. ! IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO ! BE FULL AND THE PARTIAL DERIVATIVES ARE ! STORED IN DFY AS ! DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) ! ELSE, THE JACOBIAN IS TAKEN AS BANDED AND ! THE PARTIAL DERIVATIVES ARE STORED ! DIAGONAL-WISE AS ! DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). ! ! IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN: ! IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE ! DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. ! IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC. ! ! MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: ! MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR ! ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. ! 0<=MLJAC= NUMBER OF NON-ZERO DIAGONALS BELOW ! THE MAIN DIAGONAL). ! ! MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON- ! ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). ! NEED NOT BE DEFINED IF MLJAC=N. ! ! DFX NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES ! THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X ! (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1; ! SUPPLY A DUMMY SUBROUTINE IN THE CASE IDFX=0 OR IFCN=0). ! OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM ! SUBROUTINE DFX(N,X,Y,FX) ! REAL*8 X,Y(N),FX(N) ! FX(1)= ... ! ! IDFX SWITCH FOR THE COMPUTATION OF THE DF/DX: ! IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE ! DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED. ! IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX. ! ! ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- ! ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - ! ! MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- ! MATRIX M. ! IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY ! MATRIX AND NEEDS NOT TO BE DEFINED; ! SUPPLY A DUMMY SUBROUTINE IN THIS CASE. ! IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM ! SUBROUTINE MAS(N,AM,LMAS) ! REAL*8 AM(LMAS,N) ! AM(1,1)= .... ! IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED ! AS FULL MATRIX LIKE ! AM(I,J) = M(I,J) ! ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED ! DIAGONAL-WISE AS ! AM(I-J+MUMAS+1,J) = M(I,J). ! ! IMAS GIVES INFORMATION ON THE MASS-MATRIX: ! IMAS=0: M IS SUPPOSED TO BE THE IDENTITY ! MATRIX, MAS IS NEVER CALLED. ! IMAS=1: MASS-MATRIX IS SUPPLIED. ! ! MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: ! MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ! ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. ! 0<=MLMAS= NUMBER OF NON-ZERO DIAGONALS BELOW ! THE MAIN DIAGONAL). ! MLMAS IS SUPPOSED TO BE <= MLJAC. ! ! MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- ! ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). ! NEED NOT BE DEFINED IF MLMAS=N. ! MUMAS IS SUPPOSED TO BE <= MUJAC. ! ! SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE ! NUMERICAL SOLUTION DURING INTEGRATION. ! IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. ! SUPPLY A DUMMY SUBROUTINE IF IOUT=0. ! IT MUST HAVE THE FORM ! SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN) ! REAL*8 X,Y(N) ! .... ! SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH ! GRID-POINT "X" (THEREBY THE INITIAL VALUE IS ! THE FIRST GRID-POINT). ! "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN ! IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM. ! ! IOUT GIVES INFORMATION ON THE SUBROUTINE SOLOUT: ! IOUT=0: SUBROUTINE IS NEVER CALLED ! IOUT=1: SUBROUTINE IS USED FOR OUTPUT ! ! WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". ! SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES. ! "LWORK" MUST BE AT LEAST ! N*(LJAC+LMAS+LE1+8)+5 ! WHERE ! LJAC=N IF MLJAC=N (FULL JACOBIAN) ! LJAC=MLJAC+MUJAC+1 IF MLJAC