1 SUBROUTINE DOP853(N,FCN,X,Y,XEND, 2 & RTOL,ATOL,ITOL, 3 & SOLOUT,IOUT, 4 & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) 5C ---------------------------------------------------------- 6C NUMERICAL SOLUTION OF A SYSTEM OF FIRST 0RDER 7C ORDINARY DIFFERENTIAL EQUATIONS Y'=F(X,Y). 8C THIS IS AN EXPLICIT RUNGE-KUTTA METHOD OF ORDER 8(5,3) 9C DUE TO DORMAND & PRINCE (WITH STEPSIZE CONTROL AND 10C DENSE OUTPUT) 11C 12C AUTHORS: E. HAIRER AND G. WANNER 13C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES 14C CH-1211 GENEVE 24, SWITZERLAND 15C E-MAIL: Ernst.Hairer@math.unige.ch 16C Gerhard.Wanner@math.unige.ch 17C 18C THIS CODE IS DESCRIBED IN: 19C E. HAIRER, S.P. NORSETT AND G. WANNER, SOLVING ORDINARY 20C DIFFERENTIAL EQUATIONS I. NONSTIFF PROBLEMS. 2ND EDITION. 21C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, 22C SPRINGER-VERLAG (1993) 23C 24C VERSION OF APRIL 25, 1996 25C (latest correction of a small bug: August 8, 2005) 26C 27C Edited (22 Feb 2009) by J.C. Travers: 28C renamed HINIT->HINIT853 to avoid name collision with dopri5 29C 30C INPUT PARAMETERS 31C ---------------- 32C N DIMENSION OF THE SYSTEM 33C 34C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE 35C VALUE OF F(X,Y): 36C SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR) 37C DOUBLE PRECISION X,Y(N),F(N) 38C F(1)=... ETC. 39C 40C X INITIAL X-VALUE 41C 42C Y(N) INITIAL VALUES FOR Y 43C 44C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) 45C 46C RTOL,ATOL RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY 47C CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. 48C ATOL SHOULD BE STRICTLY POSITIVE (POSSIBLY VERY SMALL) 49C 50C ITOL SWITCH FOR RTOL AND ATOL: 51C ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. 52C THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF 53C Y(I) BELOW RTOL*ABS(Y(I))+ATOL 54C ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. 55C THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW 56C RTOL(I)*ABS(Y(I))+ATOL(I). 57C 58C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE 59C NUMERICAL SOLUTION DURING INTEGRATION. 60C IF IOUT.GE.1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. 61C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. 62C IT MUST HAVE THE FORM 63C SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,CON,ICOMP,ND, 64C RPAR,IPAR,IRTRN) 65C DIMENSION Y(N),CON(8*ND),ICOMP(ND) 66C .... 67C SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH 68C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS 69C THE FIRST GRID-POINT). 70C "XOLD" IS THE PRECEDING GRID-POINT. 71C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN 72C IS SET <0, DOP853 WILL RETURN TO THE CALLING PROGRAM. 73C IF THE NUMERICAL SOLUTION IS ALTERED IN SOLOUT, 74C SET IRTRN = 2 75C 76C ----- CONTINUOUS OUTPUT: ----- 77C DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION 78C FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH 79C THE FUNCTION 80C >>> CONTD8(I,S,CON,ICOMP,ND) <<< 81C WHICH PROVIDES AN APPROXIMATION TO THE I-TH 82C COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE 83C S SHOULD LIE IN THE INTERVAL [XOLD,X]. 84C 85C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLOUT: 86C IOUT=0: SUBROUTINE IS NEVER CALLED 87C IOUT=1: SUBROUTINE IS USED FOR OUTPUT 88C IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT 89C (IN THIS CASE WORK(5) MUST BE SPECIFIED) 90C 91C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". 92C WORK(1),...,WORK(20) SERVE AS PARAMETERS FOR THE CODE. 93C FOR STANDARD USE, SET THEM TO ZERO BEFORE CALLING. 94C "LWORK" MUST BE AT LEAST 11*N+8*NRDENS+21 95C WHERE NRDENS = IWORK(5) 96C 97C LWORK DECLARED LENGTH OF ARRAY "WORK". 98C 99C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". 100C IWORK(1),...,IWORK(20) SERVE AS PARAMETERS FOR THE CODE. 101C FOR STANDARD USE, SET THEM TO ZERO BEFORE CALLING. 102C "LIWORK" MUST BE AT LEAST NRDENS+21 . 103C 104C LIWORK DECLARED LENGTH OF ARRAY "IWORK". 105C 106C RPAR, IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH 107C CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING 108C PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. 109C 110C----------------------------------------------------------------------- 111C 112C SOPHISTICATED SETTING OF PARAMETERS 113C ----------------------------------- 114C SEVERAL PARAMETERS (WORK(1),...,IWORK(1),...) ALLOW 115C TO ADAPT THE CODE TO THE PROBLEM AND TO THE NEEDS OF 116C THE USER. FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES. 117C 118C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 2.3D-16. 119C 120C WORK(2) THE SAFETY FACTOR IN STEP SIZE PREDICTION, 121C DEFAULT 0.9D0. 122C 123C WORK(3), WORK(4) PARAMETERS FOR STEP SIZE SELECTION 124C THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION 125C WORK(3) <= HNEW/HOLD <= WORK(4) 126C DEFAULT VALUES: WORK(3)=0.333D0, WORK(4)=6.D0 127C 128C WORK(5) IS THE "BETA" FOR STABILIZED STEP SIZE CONTROL 129C (SEE SECTION IV.2). POSITIVE VALUES OF BETA ( <= 0.04 ) 130C MAKE THE STEP SIZE CONTROL MORE STABLE. 131C NEGATIVE WORK(5) PROVOKE BETA=0. 132C DEFAULT 0.0D0. 133C 134C WORK(6) MAXIMAL STEP SIZE, DEFAULT XEND-X. 135C 136C WORK(7) INITIAL STEP SIZE, FOR WORK(7)=0.D0 AN INITIAL GUESS 137C IS COMPUTED WITH HELP OF THE FUNCTION HINIT 138C 139C IWORK(1) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. 140C THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000. 141C 142C IWORK(2) SWITCH FOR THE CHOICE OF THE COEFFICIENTS 143C IF IWORK(2).EQ.1 METHOD DOP853 OF DORMAND AND PRINCE 144C (SECTION II.6). 145C THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=1. 146C 147C IWORK(3) SWITCH FOR PRINTING ERROR MESSAGES 148C IF IWORK(3).LT.0 NO MESSAGES ARE BEING PRINTED 149C IF IWORK(3).GT.0 MESSAGES ARE PRINTED WITH 150C WRITE (IWORK(3),*) ... 151C DEFAULT VALUE (FOR IWORK(3)=0) IS IWORK(3)=6 152C 153C IWORK(4) TEST FOR STIFFNESS IS ACTIVATED AFTER STEP NUMBER 154C J*IWORK(4) (J INTEGER), PROVIDED IWORK(4).GT.0. 155C FOR NEGATIVE IWORK(4) THE STIFFNESS TEST IS 156C NEVER ACTIVATED; DEFAULT VALUE IS IWORK(4)=1000 157C 158C IWORK(5) = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT 159C IS REQUIRED; DEFAULT VALUE IS IWORK(5)=0; 160C FOR 0 < NRDENS < N THE COMPONENTS (FOR WHICH DENSE 161C OUTPUT IS REQUIRED) HAVE TO BE SPECIFIED IN 162C IWORK(21),...,IWORK(NRDENS+20); 163C FOR NRDENS=N THIS IS DONE BY THE CODE. 164C 165C---------------------------------------------------------------------- 166C 167C OUTPUT PARAMETERS 168C ----------------- 169C X X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED 170C (AFTER SUCCESSFUL RETURN X=XEND). 171C 172C Y(N) NUMERICAL SOLUTION AT X 173C 174C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP 175C 176C IDID REPORTS ON SUCCESSFULNESS UPON RETURN: 177C IDID= 1 COMPUTATION SUCCESSFUL, 178C IDID= 2 COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT) 179C IDID=-1 INPUT IS NOT CONSISTENT, 180C IDID=-2 LARGER NMAX IS NEEDED, 181C IDID=-3 STEP SIZE BECOMES TOO SMALL. 182C IDID=-4 PROBLEM IS PROBABLY STIFF (INTERRUPTED). 183C 184C IWORK(17) NFCN NUMBER OF FUNCTION EVALUATIONS 185C IWORK(18) NSTEP NUMBER OF COMPUTED STEPS 186C IWORK(19) NACCPT NUMBER OF ACCEPTED STEPS 187C IWORK(20) NREJCT NUMBER OF REJECTED STEPS (DUE TO ERROR TEST), 188C (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED) 189C----------------------------------------------------------------------- 190C *** *** *** *** *** *** *** *** *** *** *** *** *** 191C DECLARATIONS 192C *** *** *** *** *** *** *** *** *** *** *** *** *** 193 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 194 DIMENSION Y(N),ATOL(*),RTOL(*),WORK(LWORK),IWORK(LIWORK) 195 DIMENSION RPAR(*),IPAR(*) 196 LOGICAL ARRET 197 EXTERNAL FCN,SOLOUT 198C *** *** *** *** *** *** *** 199C SETTING THE PARAMETERS 200C *** *** *** *** *** *** *** 201 NFCN=0 202 NSTEP=0 203 NACCPT=0 204 NREJCT=0 205 ARRET=.FALSE. 206C -------- IPRINT FOR MONITORING THE PRINTING 207 IF(IWORK(3).EQ.0)THEN 208 IPRINT=6 209 ELSE 210 IPRINT=IWORK(3) 211 END IF 212C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- 213 IF(IWORK(1).EQ.0)THEN 214 NMAX=100000 215 ELSE 216 NMAX=IWORK(1) 217 IF(NMAX.LE.0)THEN 218 IF (IPRINT.GT.0) WRITE(IPRINT,*) 219 & ' WRONG INPUT IWORK(1)=',IWORK(1) 220 ARRET=.TRUE. 221 END IF 222 END IF 223C -------- METH COEFFICIENTS OF THE METHOD 224 IF(IWORK(2).EQ.0)THEN 225 METH=1 226 ELSE 227 METH=IWORK(2) 228 IF(METH.LE.0.OR.METH.GE.4)THEN 229 IF (IPRINT.GT.0) WRITE(IPRINT,*) 230 & ' CURIOUS INPUT IWORK(2)=',IWORK(2) 231 ARRET=.TRUE. 232 END IF 233 END IF 234C -------- NSTIFF PARAMETER FOR STIFFNESS DETECTION 235 NSTIFF=IWORK(4) 236 IF (NSTIFF.EQ.0) NSTIFF=1000 237 IF (NSTIFF.LT.0) NSTIFF=NMAX+10 238C -------- NRDENS NUMBER OF DENSE OUTPUT COMPONENTS 239 NRDENS=IWORK(5) 240 IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN 241 IF (IPRINT.GT.0) WRITE(IPRINT,*) 242 & ' CURIOUS INPUT IWORK(5)=',IWORK(5) 243 ARRET=.TRUE. 244 ELSE 245 IF(NRDENS.GT.0.AND.IOUT.LT.2)THEN 246 IF (IPRINT.GT.0) WRITE(IPRINT,*) 247 & ' WARNING: PUT IOUT=2 FOR DENSE OUTPUT ' 248 END IF 249 IF (NRDENS.EQ.N) THEN 250 DO I=1,NRDENS 251 IWORK(I+20)=I 252 END DO 253 END IF 254 END IF 255C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 256 IF(WORK(1).EQ.0.D0)THEN 257 UROUND=2.3D-16 258 ELSE 259 UROUND=WORK(1) 260 IF(UROUND.LE.1.D-35.OR.UROUND.GE.1.D0)THEN 261 IF (IPRINT.GT.0) WRITE(IPRINT,*) 262 & ' WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:',WORK(1) 263 ARRET=.TRUE. 264 END IF 265 END IF 266C ------- SAFETY FACTOR ------------- 267 IF(WORK(2).EQ.0.D0)THEN 268 SAFE=0.9D0 269 ELSE 270 SAFE=WORK(2) 271 IF(SAFE.GE.1.D0.OR.SAFE.LE.1.D-4)THEN 272 IF (IPRINT.GT.0) WRITE(IPRINT,*) 273 & ' CURIOUS INPUT FOR SAFETY FACTOR WORK(2)=',WORK(2) 274 ARRET=.TRUE. 275 END IF 276 END IF 277C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION 278 IF(WORK(3).EQ.0.D0)THEN 279 FAC1=0.333D0 280 ELSE 281 FAC1=WORK(3) 282 END IF 283 IF(WORK(4).EQ.0.D0)THEN 284 FAC2=6.D0 285 ELSE 286 FAC2=WORK(4) 287 END IF 288C --------- BETA FOR STEP CONTROL STABILIZATION ----------- 289 IF(WORK(5).EQ.0.D0)THEN 290 BETA=0.0D0 291 ELSE 292 IF(WORK(5).LT.0.D0)THEN 293 BETA=0.D0 294 ELSE 295 BETA=WORK(5) 296 IF(BETA.GT.0.2D0)THEN 297 IF (IPRINT.GT.0) WRITE(IPRINT,*) 298 & ' CURIOUS INPUT FOR BETA: WORK(5)=',WORK(5) 299 ARRET=.TRUE. 300 END IF 301 END IF 302 END IF 303C -------- MAXIMAL STEP SIZE 304 IF(WORK(6).EQ.0.D0)THEN 305 HMAX=XEND-X 306 ELSE 307 HMAX=WORK(6) 308 END IF 309C -------- INITIAL STEP SIZE 310 H=WORK(7) 311C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- 312 IEK1=21 313 IEK2=IEK1+N 314 IEK3=IEK2+N 315 IEK4=IEK3+N 316 IEK5=IEK4+N 317 IEK6=IEK5+N 318 IEK7=IEK6+N 319 IEK8=IEK7+N 320 IEK9=IEK8+N 321 IEK10=IEK9+N 322 IEY1=IEK10+N 323 IECO=IEY1+N 324C ------ TOTAL STORAGE REQUIREMENT ----------- 325 ISTORE=IECO+8*NRDENS-1 326 IF(ISTORE.GT.LWORK)THEN 327 IF (IPRINT.GT.0) WRITE(IPRINT,*) 328 & ' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE 329 ARRET=.TRUE. 330 END IF 331 ICOMP=21 332 ISTORE=ICOMP+NRDENS-1 333 IF(ISTORE.GT.LIWORK)THEN 334 IF (IPRINT.GT.0) WRITE(IPRINT,*) 335 & ' INSUFFICIENT STORAGE FOR IWORK, MIN. LIWORK=',ISTORE 336 ARRET=.TRUE. 337 END IF 338C -------- WHEN A FAIL HAS OCCURRED, WE RETURN WITH IDID=-1 339 IF (ARRET) THEN 340 IDID=-1 341 RETURN 342 END IF 343C -------- CALL TO CORE INTEGRATOR ------------ 344 CALL DP86CO(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IPRINT, 345 & SOLOUT,IOUT,IDID,NMAX,UROUND,METH,NSTIFF,SAFE,BETA,FAC1,FAC2, 346 & WORK(IEK1),WORK(IEK2),WORK(IEK3),WORK(IEK4),WORK(IEK5), 347 & WORK(IEK6),WORK(IEK7),WORK(IEK8),WORK(IEK9),WORK(IEK10), 348 & WORK(IEY1),WORK(IECO),IWORK(ICOMP),NRDENS,RPAR,IPAR, 349 & NFCN,NSTEP,NACCPT,NREJCT) 350 WORK(7)=H 351 IWORK(17)=NFCN 352 IWORK(18)=NSTEP 353 IWORK(19)=NACCPT 354 IWORK(20)=NREJCT 355C ----------- RETURN ----------- 356 RETURN 357 END 358C 359C 360C 361C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- 362C 363 SUBROUTINE DP86CO(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IPRINT, 364 & SOLOUT,IOUT,IDID,NMAX,UROUND,METH,NSTIFF,SAFE,BETA,FAC1,FAC2, 365 & K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,Y1,CONT,ICOMP,NRD,RPAR,IPAR, 366 & NFCN,NSTEP,NACCPT,NREJCT) 367C ---------------------------------------------------------- 368C CORE INTEGRATOR FOR DOP853 369C PARAMETERS SAME AS IN DOP853 WITH WORKSPACE ADDED 370C ---------------------------------------------------------- 371C DECLARATIONS 372C ---------------------------------------------------------- 373 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 374 parameter ( 375 & c2 = 0.526001519587677318785587544488D-01, 376 & c3 = 0.789002279381515978178381316732D-01, 377 & c4 = 0.118350341907227396726757197510D+00, 378 & c5 = 0.281649658092772603273242802490D+00, 379 & c6 = 0.333333333333333333333333333333D+00, 380 & c7 = 0.25D+00, 381 & c8 = 0.307692307692307692307692307692D+00, 382 & c9 = 0.651282051282051282051282051282D+00, 383 & c10 = 0.6D+00, 384 & c11 = 0.857142857142857142857142857142D+00, 385 & c14 = 0.1D+00, 386 & c15 = 0.2D+00, 387 & c16 = 0.777777777777777777777777777778D+00) 388 parameter ( 389 & b1 = 5.42937341165687622380535766363D-2, 390 & b6 = 4.45031289275240888144113950566D0, 391 & b7 = 1.89151789931450038304281599044D0, 392 & b8 = -5.8012039600105847814672114227D0, 393 & b9 = 3.1116436695781989440891606237D-1, 394 & b10 = -1.52160949662516078556178806805D-1, 395 & b11 = 2.01365400804030348374776537501D-1, 396 & b12 = 4.47106157277725905176885569043D-2) 397 parameter ( 398 & bhh1 = 0.244094488188976377952755905512D+00, 399 & bhh2 = 0.733846688281611857341361741547D+00, 400 & bhh3 = 0.220588235294117647058823529412D-01) 401 parameter ( 402 & er 1 = 0.1312004499419488073250102996D-01, 403 & er 6 = -0.1225156446376204440720569753D+01, 404 & er 7 = -0.4957589496572501915214079952D+00, 405 & er 8 = 0.1664377182454986536961530415D+01, 406 & er 9 = -0.3503288487499736816886487290D+00, 407 & er10 = 0.3341791187130174790297318841D+00, 408 & er11 = 0.8192320648511571246570742613D-01, 409 & er12 = -0.2235530786388629525884427845D-01) 410 parameter ( 411 & a21 = 5.26001519587677318785587544488D-2, 412 & a31 = 1.97250569845378994544595329183D-2, 413 & a32 = 5.91751709536136983633785987549D-2, 414 & a41 = 2.95875854768068491816892993775D-2, 415 & a43 = 8.87627564304205475450678981324D-2, 416 & a51 = 2.41365134159266685502369798665D-1, 417 & a53 = -8.84549479328286085344864962717D-1, 418 & a54 = 9.24834003261792003115737966543D-1, 419 & a61 = 3.7037037037037037037037037037D-2, 420 & a64 = 1.70828608729473871279604482173D-1, 421 & a65 = 1.25467687566822425016691814123D-1, 422 & a71 = 3.7109375D-2, 423 & a74 = 1.70252211019544039314978060272D-1, 424 & a75 = 6.02165389804559606850219397283D-2, 425 & a76 = -1.7578125D-2) 426 parameter ( 427 & a81 = 3.70920001185047927108779319836D-2, 428 & a84 = 1.70383925712239993810214054705D-1, 429 & a85 = 1.07262030446373284651809199168D-1, 430 & a86 = -1.53194377486244017527936158236D-2, 431 & a87 = 8.27378916381402288758473766002D-3, 432 & a91 = 6.24110958716075717114429577812D-1, 433 & a94 = -3.36089262944694129406857109825D0, 434 & a95 = -8.68219346841726006818189891453D-1, 435 & a96 = 2.75920996994467083049415600797D1, 436 & a97 = 2.01540675504778934086186788979D1, 437 & a98 = -4.34898841810699588477366255144D1, 438 & a101 = 4.77662536438264365890433908527D-1, 439 & a104 = -2.48811461997166764192642586468D0, 440 & a105 = -5.90290826836842996371446475743D-1, 441 & a106 = 2.12300514481811942347288949897D1, 442 & a107 = 1.52792336328824235832596922938D1, 443 & a108 = -3.32882109689848629194453265587D1, 444 & a109 = -2.03312017085086261358222928593D-2) 445 parameter ( 446 & a111 = -9.3714243008598732571704021658D-1, 447 & a114 = 5.18637242884406370830023853209D0, 448 & a115 = 1.09143734899672957818500254654D0, 449 & a116 = -8.14978701074692612513997267357D0, 450 & a117 = -1.85200656599969598641566180701D1, 451 & a118 = 2.27394870993505042818970056734D1, 452 & a119 = 2.49360555267965238987089396762D0, 453 & a1110 = -3.0467644718982195003823669022D0, 454 & a121 = 2.27331014751653820792359768449D0, 455 & a124 = -1.05344954667372501984066689879D1, 456 & a125 = -2.00087205822486249909675718444D0, 457 & a126 = -1.79589318631187989172765950534D1, 458 & a127 = 2.79488845294199600508499808837D1, 459 & a128 = -2.85899827713502369474065508674D0, 460 & a129 = -8.87285693353062954433549289258D0, 461 & a1210 = 1.23605671757943030647266201528D1, 462 & a1211 = 6.43392746015763530355970484046D-1) 463 parameter ( 464 & a141 = 5.61675022830479523392909219681D-2, 465 & a147 = 2.53500210216624811088794765333D-1, 466 & a148 = -2.46239037470802489917441475441D-1, 467 & a149 = -1.24191423263816360469010140626D-1, 468 & a1410 = 1.5329179827876569731206322685D-1, 469 & a1411 = 8.20105229563468988491666602057D-3, 470 & a1412 = 7.56789766054569976138603589584D-3, 471 & a1413 = -8.298D-3) 472 parameter ( 473 & a151 = 3.18346481635021405060768473261D-2, 474 & a156 = 2.83009096723667755288322961402D-2, 475 & a157 = 5.35419883074385676223797384372D-2, 476 & a158 = -5.49237485713909884646569340306D-2, 477 & a1511 = -1.08347328697249322858509316994D-4, 478 & a1512 = 3.82571090835658412954920192323D-4, 479 & a1513 = -3.40465008687404560802977114492D-4, 480 & a1514 = 1.41312443674632500278074618366D-1, 481 & a161 = -4.28896301583791923408573538692D-1, 482 & a166 = -4.69762141536116384314449447206D0, 483 & a167 = 7.68342119606259904184240953878D0, 484 & a168 = 4.06898981839711007970213554331D0, 485 & a169 = 3.56727187455281109270669543021D-1, 486 & a1613 = -1.39902416515901462129418009734D-3, 487 & a1614 = 2.9475147891527723389556272149D0, 488 & a1615 = -9.15095847217987001081870187138D0) 489 parameter ( 490 & d41 = -0.84289382761090128651353491142D+01, 491 & d46 = 0.56671495351937776962531783590D+00, 492 & d47 = -0.30689499459498916912797304727D+01, 493 & d48 = 0.23846676565120698287728149680D+01, 494 & d49 = 0.21170345824450282767155149946D+01, 495 & d410 = -0.87139158377797299206789907490D+00, 496 & d411 = 0.22404374302607882758541771650D+01, 497 & d412 = 0.63157877876946881815570249290D+00, 498 & d413 = -0.88990336451333310820698117400D-01, 499 & d414 = 0.18148505520854727256656404962D+02, 500 & d415 = -0.91946323924783554000451984436D+01, 501 & d416 = -0.44360363875948939664310572000D+01) 502 parameter ( 503 & d51 = 0.10427508642579134603413151009D+02, 504 & d56 = 0.24228349177525818288430175319D+03, 505 & d57 = 0.16520045171727028198505394887D+03, 506 & d58 = -0.37454675472269020279518312152D+03, 507 & d59 = -0.22113666853125306036270938578D+02, 508 & d510 = 0.77334326684722638389603898808D+01, 509 & d511 = -0.30674084731089398182061213626D+02, 510 & d512 = -0.93321305264302278729567221706D+01, 511 & d513 = 0.15697238121770843886131091075D+02, 512 & d514 = -0.31139403219565177677282850411D+02, 513 & d515 = -0.93529243588444783865713862664D+01, 514 & d516 = 0.35816841486394083752465898540D+02) 515 parameter ( 516 & d61 = 0.19985053242002433820987653617D+02, 517 & d66 = -0.38703730874935176555105901742D+03, 518 & d67 = -0.18917813819516756882830838328D+03, 519 & d68 = 0.52780815920542364900561016686D+03, 520 & d69 = -0.11573902539959630126141871134D+02, 521 & d610 = 0.68812326946963000169666922661D+01, 522 & d611 = -0.10006050966910838403183860980D+01, 523 & d612 = 0.77771377980534432092869265740D+00, 524 & d613 = -0.27782057523535084065932004339D+01, 525 & d614 = -0.60196695231264120758267380846D+02, 526 & d615 = 0.84320405506677161018159903784D+02, 527 & d616 = 0.11992291136182789328035130030D+02) 528 parameter ( 529 & d71 = -0.25693933462703749003312586129D+02, 530 & d76 = -0.15418974869023643374053993627D+03, 531 & d77 = -0.23152937917604549567536039109D+03, 532 & d78 = 0.35763911791061412378285349910D+03, 533 & d79 = 0.93405324183624310003907691704D+02, 534 & d710 = -0.37458323136451633156875139351D+02, 535 & d711 = 0.10409964950896230045147246184D+03, 536 & d712 = 0.29840293426660503123344363579D+02, 537 & d713 = -0.43533456590011143754432175058D+02, 538 & d714 = 0.96324553959188282948394950600D+02, 539 & d715 = -0.39177261675615439165231486172D+02, 540 & d716 = -0.14972683625798562581422125276D+03) 541 DOUBLE PRECISION Y(N),Y1(N),K1(N),K2(N),K3(N),K4(N),K5(N),K6(N) 542 DOUBLE PRECISION K7(N),K8(N),K9(N),K10(N),ATOL(*),RTOL(*) 543 DIMENSION CONT(8*NRD),ICOMP(NRD),RPAR(*),IPAR(*) 544 LOGICAL REJECT,LAST 545 EXTERNAL FCN 546 COMMON /CONDO8/XOLD,HOUT 547C *** *** *** *** *** *** *** 548C INITIALISATIONS 549C *** *** *** *** *** *** *** 550 FACOLD=1.D-4 551 EXPO1=1.d0/8.d0-BETA*0.2D0 552 FACC1=1.D0/FAC1 553 FACC2=1.D0/FAC2 554 POSNEG=SIGN(1.D0,XEND-X) 555C --- INITIAL PREPARATIONS 556 ATOLI=ATOL(1) 557 RTOLI=RTOL(1) 558 LAST=.FALSE. 559 HLAMB=0.D0 560 IASTI=0 561 CALL FCN(N,X,Y,K1,RPAR,IPAR) 562 HMAX=ABS(HMAX) 563 IORD=8 564 IF (H.EQ.0.D0) H=HINIT853(N,FCN,X,Y,XEND,POSNEG,K1,K2,K3,IORD, 565 & HMAX,ATOL,RTOL,ITOL,RPAR,IPAR) 566 NFCN=NFCN+2 567 REJECT=.FALSE. 568 XOLD=X 569 IF (IOUT.GE.1) THEN 570 IRTRN=1 571 HOUT=1.D0 572 CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,CONT,ICOMP,NRD, 573 & RPAR,IPAR,IRTRN) 574 IF (IRTRN.LT.0) GOTO 79 575 END IF 576C --- BASIC INTEGRATION STEP 577 1 CONTINUE 578 IF (NSTEP.GT.NMAX) GOTO 78 579 IF (0.1D0*ABS(H).LE.ABS(X)*UROUND)GOTO 77 580 IF ((X+1.01D0*H-XEND)*POSNEG.GT.0.D0) THEN 581 H=XEND-X 582 LAST=.TRUE. 583 END IF 584 NSTEP=NSTEP+1 585C --- THE TWELVE STAGES 586 IF (IRTRN.GE.2) THEN 587 CALL FCN(N,X,Y,K1,RPAR,IPAR) 588 END IF 589 DO 22 I=1,N 590 22 Y1(I)=Y(I)+H*A21*K1(I) 591 CALL FCN(N,X+C2*H,Y1,K2,RPAR,IPAR) 592 DO 23 I=1,N 593 23 Y1(I)=Y(I)+H*(A31*K1(I)+A32*K2(I)) 594 CALL FCN(N,X+C3*H,Y1,K3,RPAR,IPAR) 595 DO 24 I=1,N 596 24 Y1(I)=Y(I)+H*(A41*K1(I)+A43*K3(I)) 597 CALL FCN(N,X+C4*H,Y1,K4,RPAR,IPAR) 598 DO 25 I=1,N 599 25 Y1(I)=Y(I)+H*(A51*K1(I)+A53*K3(I)+A54*K4(I)) 600 CALL FCN(N,X+C5*H,Y1,K5,RPAR,IPAR) 601 DO 26 I=1,N 602 26 Y1(I)=Y(I)+H*(A61*K1(I)+A64*K4(I)+A65*K5(I)) 603 CALL FCN(N,X+C6*H,Y1,K6,RPAR,IPAR) 604 DO 27 I=1,N 605 27 Y1(I)=Y(I)+H*(A71*K1(I)+A74*K4(I)+A75*K5(I)+A76*K6(I)) 606 CALL FCN(N,X+C7*H,Y1,K7,RPAR,IPAR) 607 DO 28 I=1,N 608 28 Y1(I)=Y(I)+H*(A81*K1(I)+A84*K4(I)+A85*K5(I)+A86*K6(I)+A87*K7(I)) 609 CALL FCN(N,X+C8*H,Y1,K8,RPAR,IPAR) 610 DO 29 I=1,N 611 29 Y1(I)=Y(I)+H*(A91*K1(I)+A94*K4(I)+A95*K5(I)+A96*K6(I)+A97*K7(I) 612 & +A98*K8(I)) 613 CALL FCN(N,X+C9*H,Y1,K9,RPAR,IPAR) 614 DO 30 I=1,N 615 30 Y1(I)=Y(I)+H*(A101*K1(I)+A104*K4(I)+A105*K5(I)+A106*K6(I) 616 & +A107*K7(I)+A108*K8(I)+A109*K9(I)) 617 CALL FCN(N,X+C10*H,Y1,K10,RPAR,IPAR) 618 DO 31 I=1,N 619 31 Y1(I)=Y(I)+H*(A111*K1(I)+A114*K4(I)+A115*K5(I)+A116*K6(I) 620 & +A117*K7(I)+A118*K8(I)+A119*K9(I)+A1110*K10(I)) 621 CALL FCN(N,X+C11*H,Y1,K2,RPAR,IPAR) 622 XPH=X+H 623 DO 32 I=1,N 624 32 Y1(I)=Y(I)+H*(A121*K1(I)+A124*K4(I)+A125*K5(I)+A126*K6(I) 625 & +A127*K7(I)+A128*K8(I)+A129*K9(I)+A1210*K10(I)+A1211*K2(I)) 626 CALL FCN(N,XPH,Y1,K3,RPAR,IPAR) 627 NFCN=NFCN+11 628 DO 35 I=1,N 629 K4(I)=B1*K1(I)+B6*K6(I)+B7*K7(I)+B8*K8(I)+B9*K9(I) 630 & +B10*K10(I)+B11*K2(I)+B12*K3(I) 631 35 K5(I)=Y(I)+H*K4(I) 632C --- ERROR ESTIMATION 633 ERR=0.D0 634 ERR2=0.D0 635 IF (ITOL.EQ.0) THEN 636 DO 41 I=1,N 637 SK=ATOLI+RTOLI*MAX(ABS(Y(I)),ABS(K5(I))) 638 ERRI=K4(I)-BHH1*K1(I)-BHH2*K9(I)-BHH3*K3(I) 639 ERR2=ERR2+(ERRI/SK)**2 640 ERRI=ER1*K1(I)+ER6*K6(I)+ER7*K7(I)+ER8*K8(I)+ER9*K9(I) 641 & +ER10*K10(I)+ER11*K2(I)+ER12*K3(I) 642 41 ERR=ERR+(ERRI/SK)**2 643 ELSE 644 DO 42 I=1,N 645 SK=ATOL(I)+RTOL(I)*MAX(ABS(Y(I)),ABS(K5(I))) 646 ERRI=K4(I)-BHH1*K1(I)-BHH2*K9(I)-BHH3*K3(I) 647 ERR2=ERR2+(ERRI/SK)**2 648 ERRI=ER1*K1(I)+ER6*K6(I)+ER7*K7(I)+ER8*K8(I)+ER9*K9(I) 649 & +ER10*K10(I)+ER11*K2(I)+ER12*K3(I) 650 42 ERR=ERR+(ERRI/SK)**2 651 END IF 652 DENO=ERR+0.01D0*ERR2 653 IF (DENO.LE.0.D0) DENO=1.D0 654 ERR=ABS(H)*ERR*SQRT(1.D0/(N*DENO)) 655C --- COMPUTATION OF HNEW 656 FAC11=ERR**EXPO1 657C --- LUND-STABILIZATION 658 FAC=FAC11/FACOLD**BETA 659C --- WE REQUIRE FAC1 <= HNEW/H <= FAC2 660 FAC=MAX(FACC2,MIN(FACC1,FAC/SAFE)) 661 HNEW=H/FAC 662 IF(ERR.LE.1.D0)THEN 663C --- STEP IS ACCEPTED 664 FACOLD=MAX(ERR,1.0D-4) 665 NACCPT=NACCPT+1 666 CALL FCN(N,XPH,K5,K4,RPAR,IPAR) 667 NFCN=NFCN+1 668C ------- STIFFNESS DETECTION 669 IF (MOD(NACCPT,NSTIFF).EQ.0.OR.IASTI.GT.0) THEN 670 STNUM=0.D0 671 STDEN=0.D0 672 DO 64 I=1,N 673 STNUM=STNUM+(K4(I)-K3(I))**2 674 STDEN=STDEN+(K5(I)-Y1(I))**2 675 64 CONTINUE 676 IF (STDEN.GT.0.D0) HLAMB=ABS(H)*SQRT(STNUM/STDEN) 677 IF (HLAMB.GT.6.1D0) THEN 678 NONSTI=0 679 IASTI=IASTI+1 680 IF (IASTI.EQ.15) THEN 681 IF (IPRINT.GT.0) WRITE (IPRINT,*) 682 & ' THE PROBLEM SEEMS TO BECOME STIFF AT X = ',X 683 IF (IPRINT.LE.0) GOTO 76 684 END IF 685 ELSE 686 NONSTI=NONSTI+1 687 IF (NONSTI.EQ.6) IASTI=0 688 END IF 689 END IF 690C ------- FINAL PREPARATION FOR DENSE OUTPUT 691 IF (IOUT.GE.2) THEN 692C ---- SAVE THE FIRST FUNCTION EVALUATIONS 693 DO 62 J=1,NRD 694 I=ICOMP(J) 695 CONT(J)=Y(I) 696 YDIFF=K5(I)-Y(I) 697 CONT(J+NRD)=YDIFF 698 BSPL=H*K1(I)-YDIFF 699 CONT(J+NRD*2)=BSPL 700 CONT(J+NRD*3)=YDIFF-H*K4(I)-BSPL 701 CONT(J+NRD*4)=D41*K1(I)+D46*K6(I)+D47*K7(I)+D48*K8(I) 702 & +D49*K9(I)+D410*K10(I)+D411*K2(I)+D412*K3(I) 703 CONT(J+NRD*5)=D51*K1(I)+D56*K6(I)+D57*K7(I)+D58*K8(I) 704 & +D59*K9(I)+D510*K10(I)+D511*K2(I)+D512*K3(I) 705 CONT(J+NRD*6)=D61*K1(I)+D66*K6(I)+D67*K7(I)+D68*K8(I) 706 & +D69*K9(I)+D610*K10(I)+D611*K2(I)+D612*K3(I) 707 CONT(J+NRD*7)=D71*K1(I)+D76*K6(I)+D77*K7(I)+D78*K8(I) 708 & +D79*K9(I)+D710*K10(I)+D711*K2(I)+D712*K3(I) 709 62 CONTINUE 710C --- THE NEXT THREE FUNCTION EVALUATIONS 711 DO 51 I=1,N 712 51 Y1(I)=Y(I)+H*(A141*K1(I)+A147*K7(I)+A148*K8(I) 713 & +A149*K9(I)+A1410*K10(I)+A1411*K2(I)+A1412*K3(I) 714 & +A1413*K4(I)) 715 CALL FCN(N,X+C14*H,Y1,K10,RPAR,IPAR) 716 DO 52 I=1,N 717 52 Y1(I)=Y(I)+H*(A151*K1(I)+A156*K6(I)+A157*K7(I) 718 & +A158*K8(I)+A1511*K2(I)+A1512*K3(I)+A1513*K4(I) 719 & +A1514*K10(I)) 720 CALL FCN(N,X+C15*H,Y1,K2,RPAR,IPAR) 721 DO 53 I=1,N 722 53 Y1(I)=Y(I)+H*(A161*K1(I)+A166*K6(I)+A167*K7(I) 723 & +A168*K8(I)+A169*K9(I)+A1613*K4(I)+A1614*K10(I) 724 & +A1615*K2(I)) 725 CALL FCN(N,X+C16*H,Y1,K3,RPAR,IPAR) 726 NFCN=NFCN+3 727C --- FINAL PREPARATION 728 DO 63 J=1,NRD 729 I=ICOMP(J) 730 CONT(J+NRD*4)=H*(CONT(J+NRD*4)+D413*K4(I)+D414*K10(I) 731 & +D415*K2(I)+D416*K3(I)) 732 CONT(J+NRD*5)=H*(CONT(J+NRD*5)+D513*K4(I)+D514*K10(I) 733 & +D515*K2(I)+D516*K3(I)) 734 CONT(J+NRD*6)=H*(CONT(J+NRD*6)+D613*K4(I)+D614*K10(I) 735 & +D615*K2(I)+D616*K3(I)) 736 CONT(J+NRD*7)=H*(CONT(J+NRD*7)+D713*K4(I)+D714*K10(I) 737 & +D715*K2(I)+D716*K3(I)) 738 63 CONTINUE 739 HOUT=H 740 END IF 741 DO 67 I=1,N 742 K1(I)=K4(I) 743 67 Y(I)=K5(I) 744 XOLD=X 745 X=XPH 746 IF (IOUT.GE.1) THEN 747 CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,CONT,ICOMP,NRD, 748 & RPAR,IPAR,IRTRN) 749 IF (IRTRN.LT.0) GOTO 79 750 END IF 751C ------- NORMAL EXIT 752 IF (LAST) THEN 753 H=HNEW 754 IDID=1 755 RETURN 756 END IF 757 IF(ABS(HNEW).GT.HMAX)HNEW=POSNEG*HMAX 758 IF(REJECT)HNEW=POSNEG*MIN(ABS(HNEW),ABS(H)) 759 REJECT=.FALSE. 760 ELSE 761C --- STEP IS REJECTED 762 HNEW=H/MIN(FACC1,FAC11/SAFE) 763 REJECT=.TRUE. 764 IF(NACCPT.GE.1)NREJCT=NREJCT+1 765 LAST=.FALSE. 766 END IF 767 H=HNEW 768 GOTO 1 769C --- FAIL EXIT 770 76 CONTINUE 771 IDID=-4 772 RETURN 773 77 CONTINUE 774 IF (IPRINT.GT.0) WRITE(IPRINT,979)X 775 IF (IPRINT.GT.0) WRITE(IPRINT,*)' STEP SIZE TOO SMALL, H=',H 776 IDID=-3 777 RETURN 778 78 CONTINUE 779 IF (IPRINT.GT.0) WRITE(IPRINT,979)X 780 IF (IPRINT.GT.0) WRITE(IPRINT,*) 781 & ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED' 782 IDID=-2 783 RETURN 784 79 CONTINUE 785 IF (IPRINT.GT.0) WRITE(IPRINT,979)X 786 979 FORMAT(' EXIT OF DOP853 AT X=',E18.4) 787 IDID=2 788 RETURN 789 END 790C 791 FUNCTION HINIT853(N,FCN,X,Y,XEND,POSNEG,F0,F1,Y1,IORD, 792 & HMAX,ATOL,RTOL,ITOL,RPAR,IPAR) 793C ---------------------------------------------------------- 794C ---- COMPUTATION OF AN INITIAL STEP SIZE GUESS 795C ---------------------------------------------------------- 796 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 797 DIMENSION Y(N),Y1(N),F0(N),F1(N),ATOL(*),RTOL(*) 798 DIMENSION RPAR(*),IPAR(*) 799C ---- COMPUTE A FIRST GUESS FOR EXPLICIT EULER AS 800C ---- H = 0.01 * NORM (Y0) / NORM (F0) 801C ---- THE INCREMENT FOR EXPLICIT EULER IS SMALL 802C ---- COMPARED TO THE SOLUTION 803 DNF=0.0D0 804 DNY=0.0D0 805 ATOLI=ATOL(1) 806 RTOLI=RTOL(1) 807 IF (ITOL.EQ.0) THEN 808 DO 10 I=1,N 809 SK=ATOLI+RTOLI*ABS(Y(I)) 810 DNF=DNF+(F0(I)/SK)**2 811 10 DNY=DNY+(Y(I)/SK)**2 812 ELSE 813 DO 11 I=1,N 814 SK=ATOL(I)+RTOL(I)*ABS(Y(I)) 815 DNF=DNF+(F0(I)/SK)**2 816 11 DNY=DNY+(Y(I)/SK)**2 817 END IF 818 IF (DNF.LE.1.D-10.OR.DNY.LE.1.D-10) THEN 819 H=1.0D-6 820 ELSE 821 H=SQRT(DNY/DNF)*0.01D0 822 END IF 823 H=MIN(H,HMAX) 824 H=SIGN(H,POSNEG) 825C ---- PERFORM AN EXPLICIT EULER STEP 826 DO 12 I=1,N 827 12 Y1(I)=Y(I)+H*F0(I) 828 CALL FCN(N,X+H,Y1,F1,RPAR,IPAR) 829C ---- ESTIMATE THE SECOND DERIVATIVE OF THE SOLUTION 830 DER2=0.0D0 831 IF (ITOL.EQ.0) THEN 832 DO 15 I=1,N 833 SK=ATOLI+RTOLI*ABS(Y(I)) 834 15 DER2=DER2+((F1(I)-F0(I))/SK)**2 835 ELSE 836 DO 16 I=1,N 837 SK=ATOL(I)+RTOL(I)*ABS(Y(I)) 838 16 DER2=DER2+((F1(I)-F0(I))/SK)**2 839 END IF 840 DER2=SQRT(DER2)/H 841C ---- STEP SIZE IS COMPUTED SUCH THAT 842C ---- H**IORD * MAX ( NORM (F0), NORM (DER2)) = 0.01 843 DER12=MAX(ABS(DER2),SQRT(DNF)) 844 IF (DER12.LE.1.D-15) THEN 845 H1=MAX(1.0D-6,ABS(H)*1.0D-3) 846 ELSE 847 H1=(0.01D0/DER12)**(1.D0/IORD) 848 END IF 849 H=MIN(100*ABS(H),H1,HMAX) 850 HINIT853=SIGN(H,POSNEG) 851 RETURN 852 END 853C 854 FUNCTION CONTD8(II,X,CON,ICOMP,ND) 855C ---------------------------------------------------------- 856C THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONNECTION 857C WITH THE OUTPUT-SUBROUTINE FOR DOP853. IT PROVIDES AN 858C APPROXIMATION TO THE II-TH COMPONENT OF THE SOLUTION AT X. 859C ---------------------------------------------------------- 860 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 861 DIMENSION CON(8*ND),ICOMP(ND) 862 COMMON /CONDO8/XOLD,H 863C ----- COMPUTE PLACE OF II-TH COMPONENT 864 I=0 865 DO 5 J=1,ND 866 IF (ICOMP(J).EQ.II) I=J 867 5 CONTINUE 868 IF (I.EQ.0) THEN 869 WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',II 870 CONTD8=-1 871 RETURN 872 END IF 873 S=(X-XOLD)/H 874 S1=1.D0-S 875 CONPAR=CON(I+ND*4)+S*(CON(I+ND*5)+S1*(CON(I+ND*6)+S*CON(I+ND*7))) 876 CONTD8=CON(I)+S*(CON(I+ND)+S1*(CON(I+ND*2)+S*(CON(I+ND*3) 877 & +S1*CONPAR))) 878 RETURN 879 END 880 881