[2696] | 1 | ! DVODE_F90 demonstration program |
---|
| 2 | |
---|
| 3 | ! Toronto Nonstiff Test Set |
---|
| 4 | |
---|
| 5 | ! Reference: |
---|
| 6 | ! ALGORITHM 648, COLLECTED ALGORITHMS FROM ACM. |
---|
| 7 | ! TRANSACTIONS ON MATHEMATICAL SOFTWARE, |
---|
| 8 | ! VOL. 13, NO. 1, P. 28 |
---|
| 9 | |
---|
| 10 | ! Caution: |
---|
| 11 | ! This version of the test set is intended for use only |
---|
| 12 | ! in conjunction with this demo program for DVODE_F90. |
---|
| 13 | ! The original test routines are altered in this version. |
---|
| 14 | ! Some subroutine arguments, COMMON variables, and local |
---|
| 15 | ! variables have been moved to the following global header. |
---|
| 16 | ! You should obtain the original test suite from netlib |
---|
| 17 | ! if you plan to use the routines for another solver. |
---|
| 18 | |
---|
| 19 | !****************************************************************** |
---|
| 20 | |
---|
| 21 | MODULE NONSTIFFSET |
---|
| 22 | |
---|
| 23 | ! Note: |
---|
| 24 | ! In this version ID is defined in the main program |
---|
| 25 | ! demostiff. The other parameters are obtained by |
---|
| 26 | ! calling IVALU with a modified argument list. |
---|
| 27 | |
---|
| 28 | IMPLICIT NONE |
---|
| 29 | INTEGER IWT, N, ID, NFCN |
---|
| 30 | DOUBLE PRECISION W |
---|
| 31 | DIMENSION W(51) |
---|
| 32 | |
---|
| 33 | CONTAINS |
---|
| 34 | |
---|
| 35 | SUBROUTINE DERIVS(NEQ,T,Y,YDOT) |
---|
| 36 | IMPLICIT NONE |
---|
| 37 | INTEGER NEQ |
---|
| 38 | DOUBLE PRECISION T, Y, YDOT |
---|
| 39 | DIMENSION Y(NEQ), YDOT(NEQ) |
---|
| 40 | |
---|
| 41 | CALL FCN(T,Y,YDOT) |
---|
| 42 | RETURN |
---|
| 43 | END SUBROUTINE DERIVS |
---|
| 44 | |
---|
| 45 | SUBROUTINE IVALU(XSTART,XEND,HBEGIN,HMAX,Y) |
---|
| 46 | |
---|
| 47 | ! ROUTINE TO PROVIDE THE INITIAL VALUES REQUIRED TO SPECIFY |
---|
| 48 | ! THE MATHEMATICAL PROBLEM AS WELL AS VARIOUS PROBLEM |
---|
| 49 | ! PARAMETERS REQUIRED BY THE TESTING PACKAGE. THE APPROPRIATE |
---|
| 50 | ! SCALING VECTOR IS ALSO INITIALISED IN CASE THIS OPTION IS |
---|
| 51 | ! SELECTED. |
---|
| 52 | |
---|
| 53 | ! PARAMETERS (OUTPUT) |
---|
| 54 | ! N - DIMENSION OF THE PROBLEM |
---|
| 55 | ! XSTART - INITIAL VALUE OF THE INDEPENDENT VARIABLE |
---|
| 56 | ! XEND - FINAL VALUE OF THE INDEPENDENT VARIABLE |
---|
| 57 | ! HBEGIN - APPROPRIATE STARTING STEPSIZE |
---|
| 58 | ! Y - VECTOR OF INITIAL CONDITIONS FOR THE DEPENDENT |
---|
| 59 | ! VARIABLES |
---|
| 60 | ! WT - VECTOR OF WEIGHTS USED TO SCALE THE PROBLEM IF |
---|
| 61 | ! THIS OPTION IS SELECTED. |
---|
| 62 | |
---|
| 63 | ! PARAMETER (INPUT) |
---|
| 64 | ! IWT - FLAG TO INDICATE IF SCALED OPTION IS SELESTED |
---|
| 65 | ! ID - FLAG IDENTIFYING WHICH EQUATION IS BEING SOLVED |
---|
| 66 | |
---|
| 67 | IMPLICIT NONE |
---|
| 68 | |
---|
| 69 | ! .. Scalar Arguments .. |
---|
| 70 | DOUBLE PRECISION HBEGIN, HMAX, XEND, XSTART |
---|
| 71 | ! INTEGER ID, IWT, N |
---|
| 72 | ! .. Array Arguments .. |
---|
| 73 | DOUBLE PRECISION Y(51) |
---|
| 74 | ! DOUBLE PRECISION W(51) |
---|
| 75 | ! .. Local Scalars .. |
---|
| 76 | DOUBLE PRECISION E, HB, HM, XE, XS |
---|
| 77 | INTEGER I, IOUT |
---|
| 78 | ! .. Data statements .. |
---|
| 79 | DATA HM, HB, XS, XE/20.0D0, 1.0D0, 0.0D0, 20.0D0/ |
---|
| 80 | |
---|
| 81 | ! .. Executable Statements .. |
---|
| 82 | |
---|
| 83 | HMAX = HM |
---|
| 84 | HBEGIN = HB |
---|
| 85 | XSTART = XS |
---|
| 86 | XEND = XE |
---|
| 87 | ! GOTO (40, 60, 80, 100, 120, 20, 20, 20, 20, 20, 140, 160, 180, & |
---|
| 88 | ! 200, 220, 20, 20, 20, 20, 20, 240, 280, 320, 360, 400, 20, 20, 20,& |
---|
| 89 | ! 20, 20, 420, 420, 420, 420, 420, 20, 20, 20, 20, 20, 540, 560, & |
---|
| 90 | ! 580, 600, 620, 20, 20, 20, 20, 20, 640, 660, 680, 700, 720) ID |
---|
| 91 | GO TO (20,30,40,50,60,10,10,10,10,10,70,80,90,100,110,10,10,10,10,10, & |
---|
| 92 | 120,130,140,150,160,10,10,10,10,10,170,170,170,170,170,10,10,10,10, & |
---|
| 93 | 10,230,240,250,260,270,10,10,10,10,10,280,290,300,310,320) ID |
---|
| 94 | 10 IOUT = 6 |
---|
| 95 | WRITE (IOUT,FMT=90000) ID |
---|
| 96 | STOP |
---|
| 97 | |
---|
| 98 | ! PROBLEM CLASS A |
---|
| 99 | |
---|
| 100 | ! 1: |
---|
| 101 | 20 CONTINUE |
---|
| 102 | ! PROBLEM A1 |
---|
| 103 | N = 1 |
---|
| 104 | W(1) = 0.100D+01 |
---|
| 105 | Y(1) = 1.0D0 |
---|
| 106 | GO TO 330 |
---|
| 107 | ! 2: |
---|
| 108 | 30 CONTINUE |
---|
| 109 | ! PROBLEM A2 |
---|
| 110 | N = 1 |
---|
| 111 | W(1) = 0.100D+01 |
---|
| 112 | Y(1) = 1.0D0 |
---|
| 113 | GO TO 330 |
---|
| 114 | ! 3: |
---|
| 115 | 40 CONTINUE |
---|
| 116 | ! PROBLEM A3 |
---|
| 117 | N = 1 |
---|
| 118 | W(1) = 0.271D+01 |
---|
| 119 | Y(1) = 1.0D0 |
---|
| 120 | GO TO 330 |
---|
| 121 | ! 4: |
---|
| 122 | 50 CONTINUE |
---|
| 123 | ! PROBLEM A4 |
---|
| 124 | N = 1 |
---|
| 125 | W(1) = 0.177D+02 |
---|
| 126 | Y(1) = 1.0D0 |
---|
| 127 | GO TO 330 |
---|
| 128 | ! 5: |
---|
| 129 | 60 CONTINUE |
---|
| 130 | ! PROBLEM A5 |
---|
| 131 | N = 1 |
---|
| 132 | W(1) = 0.620D+01 |
---|
| 133 | Y(1) = 4.0D0 |
---|
| 134 | GO TO 330 |
---|
| 135 | |
---|
| 136 | ! PROBLEM CLASS B |
---|
| 137 | |
---|
| 138 | ! 11: |
---|
| 139 | 70 CONTINUE |
---|
| 140 | ! PROBLEM B1 |
---|
| 141 | N = 2 |
---|
| 142 | W(1) = 0.425D+01 |
---|
| 143 | W(2) = 0.300D+01 |
---|
| 144 | Y(1) = 1.0D0 |
---|
| 145 | Y(2) = 3.0D0 |
---|
| 146 | GO TO 330 |
---|
| 147 | ! 12: |
---|
| 148 | 80 CONTINUE |
---|
| 149 | ! PROBLEM B2 |
---|
| 150 | N = 3 |
---|
| 151 | W(1) = 0.200D+01 |
---|
| 152 | W(2) = 0.100D+01 |
---|
| 153 | W(3) = 0.100D+01 |
---|
| 154 | Y(1) = 2.0D0 |
---|
| 155 | Y(2) = 0.0D0 |
---|
| 156 | Y(3) = 1.0D0 |
---|
| 157 | GO TO 330 |
---|
| 158 | ! 13: |
---|
| 159 | 90 CONTINUE |
---|
| 160 | ! PROBLEM B3 |
---|
| 161 | N = 3 |
---|
| 162 | W(1) = 0.100D+01 |
---|
| 163 | W(2) = 0.519D+00 |
---|
| 164 | W(3) = 0.947D+00 |
---|
| 165 | Y(1) = 1.0D0 |
---|
| 166 | Y(2) = 0.0D0 |
---|
| 167 | Y(3) = 0.0D0 |
---|
| 168 | GO TO 330 |
---|
| 169 | ! 14: |
---|
| 170 | 100 CONTINUE |
---|
| 171 | ! PROBLEM B4 |
---|
| 172 | N = 3 |
---|
| 173 | W(1) = 0.300D+01 |
---|
| 174 | W(2) = 0.220D+01 |
---|
| 175 | W(3) = 0.100D+01 |
---|
| 176 | Y(1) = 3.0D0 |
---|
| 177 | Y(2) = 0.0D0 |
---|
| 178 | Y(3) = 0.0D0 |
---|
| 179 | GO TO 330 |
---|
| 180 | ! 15: |
---|
| 181 | 110 CONTINUE |
---|
| 182 | ! PROBLEM B5 |
---|
| 183 | N = 3 |
---|
| 184 | W(1) = 0.100D+01 |
---|
| 185 | W(2) = 0.100D+01 |
---|
| 186 | W(3) = 0.100D+01 |
---|
| 187 | Y(1) = 0.0D0 |
---|
| 188 | Y(2) = 1.0D0 |
---|
| 189 | Y(3) = 1.0D0 |
---|
| 190 | GO TO 330 |
---|
| 191 | |
---|
| 192 | ! PROBLEM CLASS C |
---|
| 193 | |
---|
| 194 | ! 21: |
---|
| 195 | 120 CONTINUE |
---|
| 196 | ! PROBLEM C1 |
---|
| 197 | N = 10 |
---|
| 198 | W(1) = 0.100D+01 |
---|
| 199 | W(2) = 0.368D+00 |
---|
| 200 | W(3) = 0.271D+00 |
---|
| 201 | W(4) = 0.224D+00 |
---|
| 202 | W(5) = 0.195D+00 |
---|
| 203 | W(6) = 0.175D+00 |
---|
| 204 | W(7) = 0.161D+00 |
---|
| 205 | W(8) = 0.149D+00 |
---|
| 206 | W(9) = 0.139D+00 |
---|
| 207 | W(10) = 0.998D+00 |
---|
| 208 | Y(1) = 1.0D0 |
---|
| 209 | DO I = 2, N |
---|
| 210 | Y(I) = 0.0D0 |
---|
| 211 | END DO |
---|
| 212 | GO TO 330 |
---|
| 213 | ! 22: |
---|
| 214 | 130 CONTINUE |
---|
| 215 | ! PROBLEM C2 |
---|
| 216 | N = 10 |
---|
| 217 | W(1) = 0.100D+01 |
---|
| 218 | W(2) = 0.250D+00 |
---|
| 219 | W(3) = 0.148D+00 |
---|
| 220 | W(4) = 0.105D+00 |
---|
| 221 | W(5) = 0.818D-01 |
---|
| 222 | W(6) = 0.669D-01 |
---|
| 223 | W(7) = 0.566D-01 |
---|
| 224 | W(8) = 0.491D-01 |
---|
| 225 | W(9) = 0.433D-01 |
---|
| 226 | W(10) = 0.100D+01 |
---|
| 227 | Y(1) = 1.0D0 |
---|
| 228 | DO I = 2, N |
---|
| 229 | Y(I) = 0.0D0 |
---|
| 230 | END DO |
---|
| 231 | GO TO 330 |
---|
| 232 | ! 23: |
---|
| 233 | 140 CONTINUE |
---|
| 234 | ! PROBLEM C3 |
---|
| 235 | N = 10 |
---|
| 236 | W(1) = 0.100D+01 |
---|
| 237 | W(2) = 0.204D+00 |
---|
| 238 | W(3) = 0.955D-01 |
---|
| 239 | W(4) = 0.553D-01 |
---|
| 240 | W(5) = 0.359D-01 |
---|
| 241 | W(6) = 0.252D-01 |
---|
| 242 | W(7) = 0.184D-01 |
---|
| 243 | W(8) = 0.133D-01 |
---|
| 244 | W(9) = 0.874D-02 |
---|
| 245 | W(10) = 0.435D-02 |
---|
| 246 | Y(1) = 1.0D0 |
---|
| 247 | DO I = 2, N |
---|
| 248 | Y(I) = 0.0D0 |
---|
| 249 | END DO |
---|
| 250 | GO TO 330 |
---|
| 251 | ! 24: |
---|
| 252 | 150 CONTINUE |
---|
| 253 | ! PROBLEM C4 |
---|
| 254 | N = 51 |
---|
| 255 | W(1) = 0.100D+01 |
---|
| 256 | W(2) = 0.204D+00 |
---|
| 257 | W(3) = 0.955D-01 |
---|
| 258 | W(4) = 0.553D-01 |
---|
| 259 | W(5) = 0.359D-01 |
---|
| 260 | W(6) = 0.252D-01 |
---|
| 261 | W(7) = 0.186D-01 |
---|
| 262 | W(8) = 0.143D-01 |
---|
| 263 | W(9) = 0.113D-01 |
---|
| 264 | W(10) = 0.918D-02 |
---|
| 265 | W(11) = 0.760D-02 |
---|
| 266 | W(12) = 0.622D-02 |
---|
| 267 | W(13) = 0.494D-02 |
---|
| 268 | W(14) = 0.380D-02 |
---|
| 269 | W(15) = 0.284D-02 |
---|
| 270 | W(16) = 0.207D-02 |
---|
| 271 | W(17) = 0.146D-02 |
---|
| 272 | W(18) = 0.101D-02 |
---|
| 273 | W(19) = 0.678D-03 |
---|
| 274 | W(20) = 0.444D-03 |
---|
| 275 | W(21) = 0.283D-03 |
---|
| 276 | W(22) = 0.177D-03 |
---|
| 277 | W(23) = 0.107D-03 |
---|
| 278 | W(24) = 0.637D-04 |
---|
| 279 | W(25) = 0.370D-04 |
---|
| 280 | W(26) = 0.210D-04 |
---|
| 281 | W(27) = 0.116D-04 |
---|
| 282 | W(28) = 0.631D-05 |
---|
| 283 | W(29) = 0.335D-05 |
---|
| 284 | W(30) = 0.174D-05 |
---|
| 285 | W(31) = 0.884D-06 |
---|
| 286 | W(32) = 0.440D-06 |
---|
| 287 | W(33) = 0.215D-06 |
---|
| 288 | W(34) = 0.103D-06 |
---|
| 289 | W(35) = 0.481D-07 |
---|
| 290 | W(36) = 0.221D-07 |
---|
| 291 | W(37) = 0.996D-08 |
---|
| 292 | W(38) = 0.440D-08 |
---|
| 293 | W(39) = 0.191D-08 |
---|
| 294 | W(40) = 0.814D-09 |
---|
| 295 | W(41) = 0.340D-09 |
---|
| 296 | W(42) = 0.140D-09 |
---|
| 297 | W(43) = 0.564D-10 |
---|
| 298 | W(44) = 0.224D-10 |
---|
| 299 | W(45) = 0.871D-11 |
---|
| 300 | W(46) = 0.334D-11 |
---|
| 301 | W(47) = 0.126D-11 |
---|
| 302 | W(48) = 0.465D-12 |
---|
| 303 | W(49) = 0.169D-12 |
---|
| 304 | W(50) = 0.600D-13 |
---|
| 305 | W(51) = 0.189D-13 |
---|
| 306 | Y(1) = 1.0D0 |
---|
| 307 | DO I = 2, N |
---|
| 308 | Y(I) = 0.0D0 |
---|
| 309 | END DO |
---|
| 310 | GO TO 330 |
---|
| 311 | ! 25: |
---|
| 312 | 160 CONTINUE |
---|
| 313 | ! PROBLEM C5 |
---|
| 314 | N = 30 |
---|
| 315 | W(1) = 0.545D+01 |
---|
| 316 | W(2) = 0.471D+01 |
---|
| 317 | W(3) = 0.203D+01 |
---|
| 318 | W(4) = 0.664D+01 |
---|
| 319 | W(5) = 0.834D+01 |
---|
| 320 | W(6) = 0.346D+01 |
---|
| 321 | W(7) = 0.113D+02 |
---|
| 322 | W(8) = 0.172D+02 |
---|
| 323 | W(9) = 0.748D+01 |
---|
| 324 | W(10) = 0.302D+02 |
---|
| 325 | W(11) = 0.411D+01 |
---|
| 326 | W(12) = 0.144D+01 |
---|
| 327 | W(13) = 0.244D+02 |
---|
| 328 | W(14) = 0.284D+02 |
---|
| 329 | W(15) = 0.154D+02 |
---|
| 330 | W(16) = 0.764D+00 |
---|
| 331 | W(17) = 0.661D+00 |
---|
| 332 | W(18) = 0.284D+00 |
---|
| 333 | W(19) = 0.588D+00 |
---|
| 334 | W(20) = 0.366D+00 |
---|
| 335 | W(21) = 0.169D+00 |
---|
| 336 | W(22) = 0.388D+00 |
---|
| 337 | W(23) = 0.190D+00 |
---|
| 338 | W(24) = 0.877D-01 |
---|
| 339 | W(25) = 0.413D-01 |
---|
| 340 | W(26) = 0.289D+00 |
---|
| 341 | W(27) = 0.119D+00 |
---|
| 342 | W(28) = 0.177D+00 |
---|
| 343 | W(29) = 0.246D+00 |
---|
| 344 | W(30) = 0.319D-01 |
---|
| 345 | Y(1) = 3.42947415189D0 |
---|
| 346 | Y(2) = 3.35386959711D0 |
---|
| 347 | Y(3) = 1.35494901715D0 |
---|
| 348 | Y(4) = 6.64145542550D0 |
---|
| 349 | Y(5) = 5.97156957878D0 |
---|
| 350 | Y(6) = 2.18231499728D0 |
---|
| 351 | Y(7) = 11.2630437207D0 |
---|
| 352 | Y(8) = 14.6952576794D0 |
---|
| 353 | Y(9) = 6.27960525067D0 |
---|
| 354 | Y(10) = -30.1552268759D0 |
---|
| 355 | Y(11) = 1.65699966404D0 |
---|
| 356 | Y(12) = 1.43785752721D0 |
---|
| 357 | Y(13) = -21.1238353380D0 |
---|
| 358 | Y(14) = 28.4465098142D0 |
---|
| 359 | Y(15) = 15.3882659679D0 |
---|
| 360 | Y(16) = -.557160570446D0 |
---|
| 361 | Y(17) = .505696783289D0 |
---|
| 362 | Y(18) = .230578543901D0 |
---|
| 363 | Y(19) = -.415570776342D0 |
---|
| 364 | Y(20) = .365682722812D0 |
---|
| 365 | Y(21) = .169143213293D0 |
---|
| 366 | Y(22) = -.325325669158D0 |
---|
| 367 | Y(23) = .189706021964D0 |
---|
| 368 | Y(24) = .0877265322780D0 |
---|
| 369 | Y(25) = -.0240476254170D0 |
---|
| 370 | Y(26) = -.287659532608D0 |
---|
| 371 | Y(27) = -.117219543175D0 |
---|
| 372 | Y(28) = -.176860753121D0 |
---|
| 373 | Y(29) = -.216393453025D0 |
---|
| 374 | Y(30) = -.0148647893090D0 |
---|
| 375 | GO TO 330 |
---|
| 376 | |
---|
| 377 | ! PROBLEM CLASS D |
---|
| 378 | |
---|
| 379 | 170 CONTINUE |
---|
| 380 | ! PROBLEM D1, D2, D3, D4, D5 |
---|
| 381 | E = .2D0*(REAL(ID)-3.D1) - .1D0 |
---|
| 382 | N = 4 |
---|
| 383 | IF (ID/=31) GO TO 180 |
---|
| 384 | W(1) = 0.110D+01 |
---|
| 385 | W(2) = 0.995D+00 |
---|
| 386 | W(3) = 0.101D+01 |
---|
| 387 | W(4) = 0.111D+01 |
---|
| 388 | 180 IF (ID/=32) GO TO 190 |
---|
| 389 | W(1) = 0.130D+01 |
---|
| 390 | W(2) = 0.954D+00 |
---|
| 391 | W(3) = 0.105D+01 |
---|
| 392 | W(4) = 0.136D+01 |
---|
| 393 | 190 IF (ID/=33) GO TO 200 |
---|
| 394 | W(1) = 0.150D+01 |
---|
| 395 | W(2) = 0.866D+00 |
---|
| 396 | W(3) = 0.115D+01 |
---|
| 397 | W(4) = 0.173D+01 |
---|
| 398 | 200 IF (ID/=34) GO TO 210 |
---|
| 399 | W(1) = 0.170D+01 |
---|
| 400 | W(2) = 0.714D+00 |
---|
| 401 | W(3) = 0.140D+01 |
---|
| 402 | W(4) = 0.238D+01 |
---|
| 403 | 210 IF (ID/=35) GO TO 220 |
---|
| 404 | W(1) = 0.190D+01 |
---|
| 405 | W(2) = 0.436D+00 |
---|
| 406 | W(3) = 0.229D+01 |
---|
| 407 | W(4) = 0.436D+01 |
---|
| 408 | 220 CONTINUE |
---|
| 409 | Y(1) = 1.0D0 - E |
---|
| 410 | Y(2) = 0.0D0 |
---|
| 411 | Y(3) = 0.0D0 |
---|
| 412 | Y(4) = SQRT((1.0D0+E)/(1.0D0-E)) |
---|
| 413 | GO TO 330 |
---|
| 414 | |
---|
| 415 | ! PROBLEM CLASS E |
---|
| 416 | |
---|
| 417 | ! 41: |
---|
| 418 | 230 CONTINUE |
---|
| 419 | ! PROBLEM E1 |
---|
| 420 | N = 2 |
---|
| 421 | W(1) = 0.679D+00 |
---|
| 422 | W(2) = 0.478D+00 |
---|
| 423 | E = .79788456080286536D0 |
---|
| 424 | Y(1) = E*.84147098480789651D0 |
---|
| 425 | Y(2) = E*.11956681346419146D0 |
---|
| 426 | GO TO 330 |
---|
| 427 | ! 42: |
---|
| 428 | 240 CONTINUE |
---|
| 429 | ! PROBLEM E2 |
---|
| 430 | N = 2 |
---|
| 431 | W(1) = 0.201D+01 |
---|
| 432 | W(2) = 0.268D+01 |
---|
| 433 | Y(1) = 2.0D0 |
---|
| 434 | Y(2) = 0.0D0 |
---|
| 435 | GO TO 330 |
---|
| 436 | ! 43: |
---|
| 437 | 250 CONTINUE |
---|
| 438 | ! PROBLEM E3 |
---|
| 439 | N = 2 |
---|
| 440 | W(1) = 0.116D+01 |
---|
| 441 | W(2) = 0.128D+01 |
---|
| 442 | Y(1) = 0.0D0 |
---|
| 443 | Y(2) = 0.0D0 |
---|
| 444 | GO TO 330 |
---|
| 445 | ! 44: |
---|
| 446 | 260 CONTINUE |
---|
| 447 | ! PROBLEM E4 |
---|
| 448 | N = 2 |
---|
| 449 | W(1) = 0.340D+02 |
---|
| 450 | W(2) = 0.277D+00 |
---|
| 451 | Y(1) = 3.D1 |
---|
| 452 | Y(2) = 0.D0 |
---|
| 453 | GO TO 330 |
---|
| 454 | ! 45: |
---|
| 455 | 270 CONTINUE |
---|
| 456 | ! PROBLEM E5 |
---|
| 457 | N = 2 |
---|
| 458 | W(1) = 0.141D+02 |
---|
| 459 | W(2) = 0.240D+01 |
---|
| 460 | Y(1) = 0.0D0 |
---|
| 461 | Y(2) = 0.0D0 |
---|
| 462 | GO TO 330 |
---|
| 463 | |
---|
| 464 | ! PROBLEM CLASS F |
---|
| 465 | |
---|
| 466 | ! 51: |
---|
| 467 | 280 CONTINUE |
---|
| 468 | ! PROBLEM F1 |
---|
| 469 | N = 2 |
---|
| 470 | W(1) = 0.129D+02 |
---|
| 471 | W(2) = 0.384D+02 |
---|
| 472 | Y(1) = 0.0D0 |
---|
| 473 | Y(2) = 0.0D0 |
---|
| 474 | HMAX = 1.0D0 |
---|
| 475 | GO TO 330 |
---|
| 476 | ! 52: |
---|
| 477 | 290 CONTINUE |
---|
| 478 | ! PROBLEM F2 |
---|
| 479 | N = 1 |
---|
| 480 | W(1) = 0.110D+03 |
---|
| 481 | Y(1) = 110.0D0 |
---|
| 482 | HMAX = 1.0D0 |
---|
| 483 | GO TO 330 |
---|
| 484 | ! 53: |
---|
| 485 | 300 CONTINUE |
---|
| 486 | ! PROBLEM F3 |
---|
| 487 | N = 2 |
---|
| 488 | W(1) = 0.131D+01 |
---|
| 489 | W(2) = 0.737D+00 |
---|
| 490 | Y(1) = 0.0D0 |
---|
| 491 | Y(2) = 0.0D0 |
---|
| 492 | HMAX = 1.0D0 |
---|
| 493 | HBEGIN = 0.9D0 |
---|
| 494 | GO TO 330 |
---|
| 495 | ! 54: |
---|
| 496 | 310 CONTINUE |
---|
| 497 | ! PROBLEM F4 |
---|
| 498 | N = 1 |
---|
| 499 | W(1) = 0.152D+01 |
---|
| 500 | Y(1) = 1.0D0 |
---|
| 501 | HMAX = 10.0D0 |
---|
| 502 | GO TO 330 |
---|
| 503 | ! 55: |
---|
| 504 | 320 CONTINUE |
---|
| 505 | ! PROBLEM F5 |
---|
| 506 | N = 1 |
---|
| 507 | W(1) = 0.100D+01 |
---|
| 508 | Y(1) = 1.0D0 |
---|
| 509 | HMAX = 20.0D0 |
---|
| 510 | 330 CONTINUE |
---|
| 511 | IF (IWT<0.) GO TO 340 |
---|
| 512 | DO I = 1, N |
---|
| 513 | Y(I) = Y(I)/W(I) |
---|
| 514 | END DO |
---|
| 515 | 340 CONTINUE |
---|
| 516 | RETURN |
---|
| 517 | |
---|
| 518 | 90000 FORMAT ('AN INVALID INTERNAL PROBLEM ID OF ',I4, & |
---|
| 519 | ' WAS FOUND BY THE IVALU ROUTINE',/ & |
---|
| 520 | ' RUN TERMINATED. CHECK THE DATA.') |
---|
| 521 | END SUBROUTINE IVALU |
---|
| 522 | |
---|
| 523 | SUBROUTINE EVALU(Y) |
---|
| 524 | |
---|
| 525 | ! ROUTINE TO PROVIDE THE 'TRUE' SOLUTION OF THE DIFFERENTIAL |
---|
| 526 | ! EQUATION EVALUATED AT THE ENDPOINT OF THE INTEGRATION. |
---|
| 527 | |
---|
| 528 | ! PARAMETER (OUTPUT) |
---|
| 529 | ! Y - THE TRUE SOLUTION VECTOR EVALUATED AT THE ENDPOINT |
---|
| 530 | |
---|
| 531 | ! PARAMETERS (INPUT) |
---|
| 532 | ! N - DIMENSION OF THE PROBLEM |
---|
| 533 | ! W - VECTOR OF WEIGHTS USED TO SCALE THE PROBLEM |
---|
| 534 | ! IF THIS OPTION IS SELECTED |
---|
| 535 | ! IWT - FLAG USED TO SIGNAL WHEN THE SCALED PROBLEM IS |
---|
| 536 | ! BEING SOLVED |
---|
| 537 | ! ID - FLAG USED TO INDICATE WHICH EQUATION IS BEING |
---|
| 538 | ! SOLVED |
---|
| 539 | |
---|
| 540 | IMPLICIT NONE |
---|
| 541 | |
---|
| 542 | ! .. Scalar Arguments .. |
---|
| 543 | ! INTEGER ID, IWT, N |
---|
| 544 | ! .. Array Arguments .. |
---|
| 545 | DOUBLE PRECISION Y(51) |
---|
| 546 | ! DOUBLE PRECISION W(51) |
---|
| 547 | ! .. Local Scalars .. |
---|
| 548 | INTEGER I |
---|
| 549 | |
---|
| 550 | ! .. Executable Statements .. |
---|
| 551 | |
---|
| 552 | ! GOTO (20, 40, 60, 80, 100, 620, 620, 620, 620, 620, 120, 140, 160,& |
---|
| 553 | ! 180, 200, 620, 620, 620, 620, 620, 220, 240, 260, 280, 300, 620, & |
---|
| 554 | ! 620, 620, 620, 620, 320, 340, 360, 380, 400, 620, 620, 620, 620, & |
---|
| 555 | ! 620, 420, 440, 460, 480, 500, 620, 620, 620, 620, 620, 520, 540, & |
---|
| 556 | ! 560, 580, 600) ID |
---|
| 557 | GO TO (10,20,30,40,50,310,310,310,310,310,60,70,80,90,100,310,310,310, & |
---|
| 558 | 310,310,110,120,130,140,150,310,310,310,310,310,160,170,180,190,200, & |
---|
| 559 | 310,310,310,310,310,210,220,230,240,250,310,310,310,310,310,260,270, & |
---|
| 560 | 280,290,300) ID |
---|
| 561 | GO TO 310 |
---|
| 562 | |
---|
| 563 | ! PROBLEM CLASS A |
---|
| 564 | |
---|
| 565 | ! 1: |
---|
| 566 | ! PROBLEM A1 |
---|
| 567 | 10 Y(1) = 2.061153353012535D-09 |
---|
| 568 | GO TO 310 |
---|
| 569 | ! 2: |
---|
| 570 | ! PROBLEM A2 |
---|
| 571 | 20 Y(1) = 2.182178902359887D-01 |
---|
| 572 | GO TO 310 |
---|
| 573 | ! 3: |
---|
| 574 | ! PROBLEM A3 |
---|
| 575 | 30 Y(1) = 2.491650271850414D+00 |
---|
| 576 | GO TO 310 |
---|
| 577 | ! 4: |
---|
| 578 | ! PROBLEM A4 |
---|
| 579 | 40 Y(1) = 1.773016648131483D+01 |
---|
| 580 | GO TO 310 |
---|
| 581 | ! 5: |
---|
| 582 | ! PROBLEM A5 |
---|
| 583 | 50 Y(1) = -7.887826688964196D-01 |
---|
| 584 | GO TO 310 |
---|
| 585 | |
---|
| 586 | ! PROBLEM CLASS B |
---|
| 587 | |
---|
| 588 | ! 11: |
---|
| 589 | ! PROBLEM B1 |
---|
| 590 | 60 Y(1) = 6.761876008576667D-01 |
---|
| 591 | Y(2) = 1.860816099640036D-01 |
---|
| 592 | GO TO 310 |
---|
| 593 | ! 12: |
---|
| 594 | ! PROBLEM B2 |
---|
| 595 | 70 Y(1) = 1.000000001030576D+00 |
---|
| 596 | Y(2) = 1.000000000000000D+00 |
---|
| 597 | Y(3) = 9.999999989694235D-01 |
---|
| 598 | GO TO 310 |
---|
| 599 | ! 13: |
---|
| 600 | ! PROBLEM B3 |
---|
| 601 | 80 Y(1) = 2.061153488557776D-09 |
---|
| 602 | Y(2) = 5.257228022048349D-02 |
---|
| 603 | Y(3) = 9.474277177183630D-01 |
---|
| 604 | GO TO 310 |
---|
| 605 | ! 14: |
---|
| 606 | ! PROBLEM B4 |
---|
| 607 | 90 Y(1) = 9.826950928005993D-01 |
---|
| 608 | Y(2) = 2.198447081694832D+00 |
---|
| 609 | Y(3) = 9.129452507276399D-01 |
---|
| 610 | GO TO 310 |
---|
| 611 | ! 15: |
---|
| 612 | ! PROBLEM B5 |
---|
| 613 | 100 Y(1) = -9.396570798729192D-01 |
---|
| 614 | Y(2) = -3.421177754000779D-01 |
---|
| 615 | Y(3) = 7.414126596199957D-01 |
---|
| 616 | GO TO 310 |
---|
| 617 | |
---|
| 618 | ! PROBLEM CLASS C |
---|
| 619 | |
---|
| 620 | ! 21: |
---|
| 621 | ! PROBLEM C1 |
---|
| 622 | 110 Y(1) = 2.061153622240064D-09 |
---|
| 623 | Y(2) = 4.122307244619555D-08 |
---|
| 624 | Y(3) = 4.122307244716968D-07 |
---|
| 625 | Y(4) = 2.748204829855288D-06 |
---|
| 626 | Y(5) = 1.374102414941961D-05 |
---|
| 627 | Y(6) = 5.496409659803266D-05 |
---|
| 628 | Y(7) = 1.832136553274552D-04 |
---|
| 629 | Y(8) = 5.234675866508716D-04 |
---|
| 630 | Y(9) = 1.308668966628220D-03 |
---|
| 631 | Y(10) = 9.979127409508656D-01 |
---|
| 632 | GO TO 310 |
---|
| 633 | ! 22: |
---|
| 634 | ! PROBLEM C2 |
---|
| 635 | 120 Y(1) = 2.061153577984930D-09 |
---|
| 636 | Y(2) = 2.061153573736588D-09 |
---|
| 637 | Y(3) = 2.061153569488245D-09 |
---|
| 638 | Y(4) = 2.061153565239902D-09 |
---|
| 639 | Y(5) = 2.061153560991560D-09 |
---|
| 640 | Y(6) = 2.061153556743217D-09 |
---|
| 641 | Y(7) = 2.061153552494874D-09 |
---|
| 642 | Y(8) = 2.061153548246532D-09 |
---|
| 643 | Y(9) = 2.061153543998189D-09 |
---|
| 644 | Y(10) = 9.999999814496180D-01 |
---|
| 645 | GO TO 310 |
---|
| 646 | ! 23: |
---|
| 647 | ! PROBLEM C3 |
---|
| 648 | 130 Y(1) = 2.948119211022058D-03 |
---|
| 649 | Y(2) = 5.635380154844266D-03 |
---|
| 650 | Y(3) = 7.829072515926013D-03 |
---|
| 651 | Y(4) = 9.348257908594937D-03 |
---|
| 652 | Y(5) = 1.007943610301970D-02 |
---|
| 653 | Y(6) = 9.982674171429909D-03 |
---|
| 654 | Y(7) = 9.088693332766085D-03 |
---|
| 655 | Y(8) = 7.489115195185912D-03 |
---|
| 656 | Y(9) = 5.322964130953349D-03 |
---|
| 657 | Y(10) = 2.762434379029886D-03 |
---|
| 658 | GO TO 310 |
---|
| 659 | ! 24: |
---|
| 660 | ! PROBLEM C4 |
---|
| 661 | 140 Y(1) = 3.124111453721466D-03 |
---|
| 662 | Y(2) = 6.015416842150318D-03 |
---|
| 663 | Y(3) = 8.470021834842650D-03 |
---|
| 664 | Y(4) = 1.033682931733337D-02 |
---|
| 665 | Y(5) = 1.153249572873923D-02 |
---|
| 666 | Y(6) = 1.204549525737964D-02 |
---|
| 667 | Y(7) = 1.192957068015293D-02 |
---|
| 668 | Y(8) = 1.128883207111195D-02 |
---|
| 669 | Y(9) = 1.025804501391024D-02 |
---|
| 670 | Y(10) = 8.982017581934167D-03 |
---|
| 671 | Y(11) = 7.597500902492453D-03 |
---|
| 672 | Y(12) = 6.219920556824985D-03 |
---|
| 673 | Y(13) = 4.935916341009131D-03 |
---|
| 674 | Y(14) = 3.801432544256119D-03 |
---|
| 675 | Y(15) = 2.844213677587894D-03 |
---|
| 676 | Y(16) = 2.069123394222672D-03 |
---|
| 677 | Y(17) = 1.464687282843915D-03 |
---|
| 678 | Y(18) = 1.009545263941126D-03 |
---|
| 679 | Y(19) = 6.779354330227017D-04 |
---|
| 680 | Y(20) = 4.437815269118510D-04 |
---|
| 681 | Y(21) = 2.833264542938954D-04 |
---|
| 682 | Y(22) = 1.765005798796805D-04 |
---|
| 683 | Y(23) = 1.073342592697238D-04 |
---|
| 684 | Y(24) = 6.374497601777217D-05 |
---|
| 685 | Y(25) = 3.698645309704183D-05 |
---|
| 686 | Y(26) = 2.097466832643746D-05 |
---|
| 687 | Y(27) = 1.162956710412555D-05 |
---|
| 688 | Y(28) = 6.306710405783322D-06 |
---|
| 689 | Y(29) = 3.346286430868515D-06 |
---|
| 690 | Y(30) = 1.737760074184334D-06 |
---|
| 691 | Y(31) = 8.835366904275847D-07 |
---|
| 692 | Y(32) = 4.399520411127637D-07 |
---|
| 693 | Y(33) = 2.146181897152360D-07 |
---|
| 694 | Y(34) = 1.025981211654928D-07 |
---|
| 695 | Y(35) = 4.807864068784215D-08 |
---|
| 696 | Y(36) = 2.209175152474847D-08 |
---|
| 697 | Y(37) = 9.956251263138180D-09 |
---|
| 698 | Y(38) = 4.402193653748924D-09 |
---|
| 699 | Y(39) = 1.910149382204028D-09 |
---|
| 700 | Y(40) = 8.135892921473050D-10 |
---|
| 701 | Y(41) = 3.402477118549235D-10 |
---|
| 702 | Y(42) = 1.397485617545782D-10 |
---|
| 703 | Y(43) = 5.638575303049199D-11 |
---|
| 704 | Y(44) = 2.235459707956947D-11 |
---|
| 705 | Y(45) = 8.710498036398032D-12 |
---|
| 706 | Y(46) = 3.336554275346643D-12 |
---|
| 707 | Y(47) = 1.256679567784939D-12 |
---|
| 708 | Y(48) = 4.654359053128788D-13 |
---|
| 709 | Y(49) = 1.693559145599857D-13 |
---|
| 710 | Y(50) = 5.996593816663054D-14 |
---|
| 711 | Y(51) = 1.891330702629865D-14 |
---|
| 712 | GO TO 310 |
---|
| 713 | ! 25: |
---|
| 714 | ! PROBLEM C5 |
---|
| 715 | 150 Y(1) = -4.792730224323733D+00 |
---|
| 716 | Y(2) = -2.420550725448973D+00 |
---|
| 717 | Y(3) = -9.212509306014886D-01 |
---|
| 718 | Y(4) = -4.217310404035213D+00 |
---|
| 719 | Y(5) = 7.356202947498970D+00 |
---|
| 720 | Y(6) = 3.223785985421212D+00 |
---|
| 721 | Y(7) = 4.035559443262270D+00 |
---|
| 722 | Y(8) = 1.719865528670555D+01 |
---|
| 723 | Y(9) = 7.478910794233703D+00 |
---|
| 724 | Y(10) = -2.998759326324844D+01 |
---|
| 725 | Y(11) = -4.107310937550929D+00 |
---|
| 726 | Y(12) = -9.277008321754407D-01 |
---|
| 727 | Y(13) = -2.442125302518482D+01 |
---|
| 728 | Y(14) = 2.381459045746554D+01 |
---|
| 729 | Y(15) = 1.492096306951359D+01 |
---|
| 730 | Y(16) = 3.499208963063806D-01 |
---|
| 731 | Y(17) = -5.748487687912825D-01 |
---|
| 732 | Y(18) = -2.551694020879149D-01 |
---|
| 733 | Y(19) = -5.237040978903326D-01 |
---|
| 734 | Y(20) = -2.493000463579661D-01 |
---|
| 735 | Y(21) = -8.045341642044464D-02 |
---|
| 736 | Y(22) = -3.875289237334110D-01 |
---|
| 737 | Y(23) = 5.648603288767891D-02 |
---|
| 738 | Y(24) = 3.023606472143342D-02 |
---|
| 739 | Y(25) = 4.133856546712445D-02 |
---|
| 740 | Y(26) = -2.862393029841379D-01 |
---|
| 741 | Y(27) = -1.183032405136207D-01 |
---|
| 742 | Y(28) = -1.511986457359206D-01 |
---|
| 743 | Y(29) = -2.460068894318766D-01 |
---|
| 744 | Y(30) = -3.189687411323877D-02 |
---|
| 745 | GO TO 310 |
---|
| 746 | |
---|
| 747 | ! PROBLEM CLASS D |
---|
| 748 | |
---|
| 749 | ! 31: |
---|
| 750 | ! PROBLEM D1 |
---|
| 751 | 160 Y(1) = 2.198835352008397D-01 |
---|
| 752 | Y(2) = 9.427076846341813D-01 |
---|
| 753 | Y(3) = -9.787659841058176D-01 |
---|
| 754 | Y(4) = 3.287977990962036D-01 |
---|
| 755 | GO TO 310 |
---|
| 756 | ! 32: |
---|
| 757 | ! PROBLEM D2 |
---|
| 758 | 170 Y(1) = -1.777027357140412D-01 |
---|
| 759 | Y(2) = 9.467784719905892D-01 |
---|
| 760 | Y(3) = -1.030294163192969D+00 |
---|
| 761 | Y(4) = 1.211074890053952D-01 |
---|
| 762 | GO TO 310 |
---|
| 763 | ! 33: |
---|
| 764 | ! PROBLEM D3 |
---|
| 765 | 180 Y(1) = -5.780432953035361D-01 |
---|
| 766 | Y(2) = 8.633840009194193D-01 |
---|
| 767 | Y(3) = -9.595083730380727D-01 |
---|
| 768 | Y(4) = -6.504915126712089D-02 |
---|
| 769 | GO TO 310 |
---|
| 770 | ! 34: |
---|
| 771 | ! PROBLEM D4 |
---|
| 772 | 190 Y(1) = -9.538990293416394D-01 |
---|
| 773 | Y(2) = 6.907409024219432D-01 |
---|
| 774 | Y(3) = -8.212674270877433D-01 |
---|
| 775 | Y(4) = -1.539574259125825D-01 |
---|
| 776 | GO TO 310 |
---|
| 777 | ! 35: |
---|
| 778 | ! PROBLEM D5 |
---|
| 779 | 200 Y(1) = -1.295266250987574D+00 |
---|
| 780 | Y(2) = 4.003938963792321D-01 |
---|
| 781 | Y(3) = -6.775390924707566D-01 |
---|
| 782 | Y(4) = -1.270838154278686D-01 |
---|
| 783 | GO TO 310 |
---|
| 784 | |
---|
| 785 | ! PROBLEM CLASS E |
---|
| 786 | |
---|
| 787 | ! 41: |
---|
| 788 | ! PROBLEM E1 |
---|
| 789 | 210 Y(1) = 1.456723600728308D-01 |
---|
| 790 | Y(2) = -9.883500195574063D-02 |
---|
| 791 | GO TO 310 |
---|
| 792 | ! 42: |
---|
| 793 | ! PROBLEM E2 |
---|
| 794 | 220 Y(1) = 2.008149762174948D+00 |
---|
| 795 | Y(2) = -4.250887527320057D-02 |
---|
| 796 | GO TO 310 |
---|
| 797 | ! 43: |
---|
| 798 | ! PROBLEM E3 |
---|
| 799 | 230 Y(1) = -1.004178858647128D-01 |
---|
| 800 | Y(2) = 2.411400132095954D-01 |
---|
| 801 | GO TO 310 |
---|
| 802 | ! 44: |
---|
| 803 | ! PROBLEM E4 |
---|
| 804 | 240 Y(1) = 3.395091444646555D+01 |
---|
| 805 | Y(2) = 2.767822659672869D-01 |
---|
| 806 | GO TO 310 |
---|
| 807 | ! 45: |
---|
| 808 | ! PROBLEM E5 |
---|
| 809 | 250 Y(1) = 1.411797390542629D+01 |
---|
| 810 | Y(2) = 2.400000000000002D+00 |
---|
| 811 | GO TO 310 |
---|
| 812 | |
---|
| 813 | ! PROBLEM CLASS F |
---|
| 814 | |
---|
| 815 | ! 51: |
---|
| 816 | ! PROBLEM F1 |
---|
| 817 | 260 Y(1) = -1.294460621213470D1 |
---|
| 818 | Y(2) = -2.208575158908672D-15 |
---|
| 819 | GO TO 310 |
---|
| 820 | ! 52: |
---|
| 821 | ! PROBLEM F2 |
---|
| 822 | 270 Y(1) = 70.03731057008607D0 |
---|
| 823 | GO TO 310 |
---|
| 824 | ! 53: |
---|
| 825 | ! PROBLEM F3 |
---|
| 826 | 280 Y(1) = -3.726957553088175D-1 |
---|
| 827 | Y(2) = -6.230137949234190D-1 |
---|
| 828 | GO TO 310 |
---|
| 829 | ! 54: |
---|
| 830 | ! PROBLEM F4 |
---|
| 831 | 290 Y(1) = 9.815017249707434D-11 |
---|
| 832 | GO TO 310 |
---|
| 833 | ! 55: |
---|
| 834 | ! PROBLEM F5 |
---|
| 835 | 300 Y(1) = 1.0D0 |
---|
| 836 | 310 CONTINUE |
---|
| 837 | IF (IWT<0) GO TO 320 |
---|
| 838 | DO I = 1, N |
---|
| 839 | Y(I) = Y(I)/W(I) |
---|
| 840 | END DO |
---|
| 841 | 320 CONTINUE |
---|
| 842 | RETURN |
---|
| 843 | END SUBROUTINE EVALU |
---|
| 844 | |
---|
| 845 | SUBROUTINE FCN(X,Y,YP) |
---|
| 846 | |
---|
| 847 | ! ROUTINE TO EVALUATE THE DERIVATIVE F(X,Y) CORRESPONDING TO |
---|
| 848 | ! THE DIFFERENTIAL EQUATION: |
---|
| 849 | ! DY/DX = F(X,Y) . |
---|
| 850 | ! THE ROUTINE STORES THE VECTOR OF DERIVATIVES IN YP(*). THE |
---|
| 851 | ! PARTICULAR EQUATION BEING INTEGRATED IS INDICATED BY THE |
---|
| 852 | ! VALUE OF THE FLAG ID WHICH IS PASSED THROUGH COMMON. THE |
---|
| 853 | ! DIFFERENTIAL EQUATION IS SCALED BY THE WEIGHT VECTOR W(*) |
---|
| 854 | ! IF THIS OPTION HAS BEEN SELECTED (IF SO IT IS SIGNALLED |
---|
| 855 | ! BY THE FLAG IWT). |
---|
| 856 | |
---|
| 857 | IMPLICIT NONE |
---|
| 858 | |
---|
| 859 | ! .. Scalar Arguments .. |
---|
| 860 | DOUBLE PRECISION X |
---|
| 861 | ! .. Array Arguments .. |
---|
| 862 | DOUBLE PRECISION Y(51), YP(51) |
---|
| 863 | ! .. Local Scalars .. |
---|
| 864 | DOUBLE PRECISION C1, C2, D, EX, K2, M0, P, TEMP |
---|
| 865 | INTEGER I, I3, I3M2, ITEMP, J, L, LL, MM |
---|
| 866 | ! .. Local Arrays .. |
---|
| 867 | DOUBLE PRECISION M(5), Q(5,5), R(5), YTEMP(51) |
---|
| 868 | ! .. Data statements .. |
---|
| 869 | ! THE FOLLOWING DATA IS FOR PROBLEM C5 AND DEFINES THE MASSES |
---|
| 870 | ! OF THE 5 OUTER PLANETS ETC. IN SOLAR UNITS. |
---|
| 871 | ! K2 IS THE GRAVITATIONAL CONSTANT. |
---|
| 872 | ! THE NEXT DATA IS FOR PROBLEMS F1 AND F5. |
---|
| 873 | ! C1 IS PI**2 + 0.1**2 AND C2 IS SUM I**(4/3) FOR I=1 TO 19. |
---|
| 874 | DATA M0/1.00000597682D0/, M/.954786104043D-3, .285583733151D-3, & |
---|
| 875 | .437273164546D-4, .517759138449D-4, .277777777778D-5/ |
---|
| 876 | DATA K2/2.95912208286D0/ |
---|
| 877 | DATA EX, C1, C2/.33333333333333333D0, 9.879604401089358D0, & |
---|
| 878 | 438.4461015267790D0/ |
---|
| 879 | |
---|
| 880 | ! .. Executable Statements .. |
---|
| 881 | |
---|
| 882 | NFCN = NFCN + 1 |
---|
| 883 | IF (IWT<0) GO TO 10 |
---|
| 884 | DO I = 1, N |
---|
| 885 | YTEMP(I) = Y(I) |
---|
| 886 | Y(I) = Y(I)*W(I) |
---|
| 887 | END DO |
---|
| 888 | 10 CONTINUE |
---|
| 889 | GO TO (20,30,40,50,60,330,330,330,330,330,70,80,90,100,110,330,330, & |
---|
| 890 | 330,330,330,120,130,140,150,160,330,330,330,330,330,190,190,190,190, & |
---|
| 891 | 190,330,330,330,330,330,200,210,220,230,240,330,330,330,330,330,250, & |
---|
| 892 | 270,290,300,320) ID |
---|
| 893 | GO TO 330 |
---|
| 894 | |
---|
| 895 | ! PROBLEM CLASS A |
---|
| 896 | |
---|
| 897 | ! 1: |
---|
| 898 | ! PROBLEM A1 |
---|
| 899 | 20 YP(1) = -Y(1) |
---|
| 900 | GO TO 330 |
---|
| 901 | ! 2: |
---|
| 902 | ! PROBLEM A2 |
---|
| 903 | 30 YP(1) = -.5D0*Y(1)*Y(1)*Y(1) |
---|
| 904 | GO TO 330 |
---|
| 905 | ! 3: |
---|
| 906 | ! PROBLEM A3 |
---|
| 907 | 40 YP(1) = Y(1)*COS(X) |
---|
| 908 | GO TO 330 |
---|
| 909 | ! 4: |
---|
| 910 | ! PROBLEM A4 |
---|
| 911 | 50 YP(1) = (1.0D0-Y(1)/20.0D0)*Y(1)/4.D0 |
---|
| 912 | GO TO 330 |
---|
| 913 | ! 5: |
---|
| 914 | ! PROBLEM A5 |
---|
| 915 | 60 YP(1) = (Y(1)-X)/(Y(1)+X) |
---|
| 916 | GO TO 330 |
---|
| 917 | |
---|
| 918 | ! PROBLEM CLASS B |
---|
| 919 | |
---|
| 920 | ! 11: |
---|
| 921 | ! PROBLEM B1 |
---|
| 922 | 70 D = Y(1) - Y(1)*Y(2) |
---|
| 923 | YP(1) = D + D |
---|
| 924 | YP(2) = -(Y(2)-Y(1)*Y(2)) |
---|
| 925 | GO TO 330 |
---|
| 926 | ! 12: |
---|
| 927 | ! PROBLEM B2 |
---|
| 928 | 80 YP(1) = -Y(1) + Y(2) |
---|
| 929 | YP(3) = Y(2) - Y(3) |
---|
| 930 | YP(2) = -YP(1) - YP(3) |
---|
| 931 | GO TO 330 |
---|
| 932 | ! 13: |
---|
| 933 | ! PROBLEM B3 |
---|
| 934 | 90 D = Y(2)*Y(2) |
---|
| 935 | YP(1) = -Y(1) |
---|
| 936 | YP(2) = Y(1) - D |
---|
| 937 | YP(3) = D |
---|
| 938 | GO TO 330 |
---|
| 939 | ! 14: |
---|
| 940 | ! PROBLEM B4 |
---|
| 941 | 100 D = SQRT(Y(1)*Y(1)+Y(2)*Y(2)) |
---|
| 942 | YP(1) = -Y(2) - Y(1)*Y(3)/D |
---|
| 943 | YP(2) = Y(1) - Y(2)*Y(3)/D |
---|
| 944 | YP(3) = Y(1)/D |
---|
| 945 | GO TO 330 |
---|
| 946 | ! 15: |
---|
| 947 | ! PROBLEM B5 |
---|
| 948 | 110 YP(1) = Y(2)*Y(3) |
---|
| 949 | YP(2) = -Y(1)*Y(3) |
---|
| 950 | YP(3) = -.51D0*Y(1)*Y(2) |
---|
| 951 | GO TO 330 |
---|
| 952 | |
---|
| 953 | ! PROBLEM CLASS C |
---|
| 954 | |
---|
| 955 | ! 21: |
---|
| 956 | ! PROBLEM C1 |
---|
| 957 | 120 YP(1) = -Y(1) |
---|
| 958 | DO I = 2, 9 |
---|
| 959 | YP(I) = Y(I-1) - Y(I) |
---|
| 960 | END DO |
---|
| 961 | YP(10) = Y(9) |
---|
| 962 | GO TO 330 |
---|
| 963 | ! 22: |
---|
| 964 | ! PROBLEM C2 |
---|
| 965 | 130 YP(1) = -Y(1) |
---|
| 966 | DO I = 2, 9 |
---|
| 967 | YP(I) = REAL(I-1)*Y(I-1) - REAL(I)*Y(I) |
---|
| 968 | END DO |
---|
| 969 | YP(10) = 9.0D0*Y(9) |
---|
| 970 | GO TO 330 |
---|
| 971 | ! 23: |
---|
| 972 | ! PROBLEM C3 |
---|
| 973 | 140 YP(1) = -2.0D0*Y(1) + Y(2) |
---|
| 974 | DO I = 2, 9 |
---|
| 975 | YP(I) = Y(I-1) - 2.0D0*Y(I) + Y(I+1) |
---|
| 976 | END DO |
---|
| 977 | YP(10) = Y(9) - 2.0D0*Y(10) |
---|
| 978 | GO TO 330 |
---|
| 979 | ! 24: |
---|
| 980 | ! PROBLEM C4 |
---|
| 981 | 150 YP(1) = -2.0D0*Y(1) + Y(2) |
---|
| 982 | DO I = 2, 50 |
---|
| 983 | YP(I) = Y(I-1) - 2.0D0*Y(I) + Y(I+1) |
---|
| 984 | END DO |
---|
| 985 | YP(51) = Y(50) - 2.0D0*Y(51) |
---|
| 986 | GO TO 330 |
---|
| 987 | ! 25: |
---|
| 988 | ! PROBLEM C5 |
---|
| 989 | 160 I = 0 |
---|
| 990 | DO L = 3, 15, 3 |
---|
| 991 | I = I + 1 |
---|
| 992 | P = Y(L-2)**2 + Y(L-1)**2 + Y(L)**2 |
---|
| 993 | R(I) = 1.0D0/(P*SQRT(P)) |
---|
| 994 | J = 0 |
---|
| 995 | DO LL = 3, 15, 3 |
---|
| 996 | J = J + 1 |
---|
| 997 | IF (LL/=L) GO TO 170 |
---|
| 998 | GO TO 180 |
---|
| 999 | ! THEN |
---|
| 1000 | 170 P = (Y(L-2)-Y(LL-2))**2 + (Y(L-1)-Y(LL-1))**2 + (Y(L)-Y(LL))**2 |
---|
| 1001 | Q(I,J) = 1.0D0/(P*SQRT(P)) |
---|
| 1002 | Q(J,I) = Q(I,J) |
---|
| 1003 | 180 CONTINUE |
---|
| 1004 | END DO |
---|
| 1005 | END DO |
---|
| 1006 | I3 = 0 |
---|
| 1007 | DO I = 1, 5 |
---|
| 1008 | I3 = I3 + 3 |
---|
| 1009 | I3M2 = I3 - 2 |
---|
| 1010 | DO LL = I3M2, I3 |
---|
| 1011 | MM = LL - I3 |
---|
| 1012 | YP(LL) = Y(LL+15) |
---|
| 1013 | P = 0.0D0 |
---|
| 1014 | DO J = 1, 5 |
---|
| 1015 | MM = MM + 3 |
---|
| 1016 | IF (J/=I) P = P + M(J)*(Y(MM)*(Q(I,J)-R(J))-Y(LL)*Q(I,J)) |
---|
| 1017 | END DO |
---|
| 1018 | YP(LL+15) = K2*(-(M0+M(I))*Y(LL)*R(I)+P) |
---|
| 1019 | END DO |
---|
| 1020 | END DO |
---|
| 1021 | GO TO 330 |
---|
| 1022 | |
---|
| 1023 | ! PROBLEM CLASS D |
---|
| 1024 | |
---|
| 1025 | ! 31:32:33:34:35: |
---|
| 1026 | ! PROBLEMS D1, D2, D3, D4, D5 |
---|
| 1027 | 190 YP(1) = Y(3) |
---|
| 1028 | YP(2) = Y(4) |
---|
| 1029 | D = Y(1)*Y(1) + Y(2)*Y(2) |
---|
| 1030 | D = SQRT(D*D*D) |
---|
| 1031 | YP(3) = -Y(1)/D |
---|
| 1032 | YP(4) = -Y(2)/D |
---|
| 1033 | GO TO 330 |
---|
| 1034 | |
---|
| 1035 | ! PROBLEM CLASS E |
---|
| 1036 | |
---|
| 1037 | ! 41: |
---|
| 1038 | ! PROBLEM E1 |
---|
| 1039 | 200 YP(1) = Y(2) |
---|
| 1040 | YP(2) = -(Y(2)/(X+1.0D0)+(1.0D0-.25D0/(X+1.0D0)**2)*Y(1)) |
---|
| 1041 | GO TO 330 |
---|
| 1042 | ! 42: |
---|
| 1043 | ! PROBLEM E2 |
---|
| 1044 | 210 YP(1) = Y(2) |
---|
| 1045 | YP(2) = (1.0D0-Y(1)*Y(1))*Y(2) - Y(1) |
---|
| 1046 | GO TO 330 |
---|
| 1047 | ! 43: |
---|
| 1048 | ! PROBLEM E3 |
---|
| 1049 | 220 YP(1) = Y(2) |
---|
| 1050 | YP(2) = Y(1)**3/6.0D0 - Y(1) + 2.0D0*SIN(2.78535D0*X) |
---|
| 1051 | GO TO 330 |
---|
| 1052 | ! 44: |
---|
| 1053 | ! PROBLEM E4 |
---|
| 1054 | 230 YP(1) = Y(2) |
---|
| 1055 | YP(2) = .032D0 - .4D0*Y(2)*Y(2) |
---|
| 1056 | GO TO 330 |
---|
| 1057 | ! 45: |
---|
| 1058 | ! PROBLEM E5 |
---|
| 1059 | 240 YP(1) = Y(2) |
---|
| 1060 | YP(2) = SQRT(1.0D0+Y(2)*Y(2))/(25.0D0-X) |
---|
| 1061 | GO TO 330 |
---|
| 1062 | |
---|
| 1063 | ! PROBLEM CLASS F |
---|
| 1064 | |
---|
| 1065 | ! 51: |
---|
| 1066 | ! PROBLEM F1 |
---|
| 1067 | 250 YP(1) = Y(2) |
---|
| 1068 | YP(2) = .2D0*Y(2) - C1*Y(1) |
---|
| 1069 | ITEMP = INT(X) |
---|
| 1070 | IF ((ITEMP/2)*2==ITEMP) GO TO 260 |
---|
| 1071 | YP(2) = YP(2) - 1.0D0 |
---|
| 1072 | GO TO 330 |
---|
| 1073 | 260 YP(2) = YP(2) + 1.0D0 |
---|
| 1074 | GO TO 330 |
---|
| 1075 | ! 52: |
---|
| 1076 | ! PROBLEM F2 |
---|
| 1077 | 270 ITEMP = INT(X) |
---|
| 1078 | IF ((ITEMP/2)*2==ITEMP) GO TO 280 |
---|
| 1079 | YP(1) = 55.0D0 - .50D0*Y(1) |
---|
| 1080 | GO TO 330 |
---|
| 1081 | 280 YP(1) = 55.0D0 - 1.50D0*Y(1) |
---|
| 1082 | GO TO 330 |
---|
| 1083 | ! 53: |
---|
| 1084 | ! PROBLEM F3 |
---|
| 1085 | 290 YP(1) = Y(2) |
---|
| 1086 | YP(2) = .01D0*Y(2)*(1.0D0-Y(1)**2) - Y(1) - ABS(SIN( & |
---|
| 1087 | 3.1415926535897932D0*X)) |
---|
| 1088 | GO TO 330 |
---|
| 1089 | ! 54: |
---|
| 1090 | ! PROBLEM F4 |
---|
| 1091 | 300 IF (X>10.) GO TO 310 |
---|
| 1092 | TEMP = X - 5.0D0 |
---|
| 1093 | YP(1) = -2.0D0/21.0D0 - 120.0D0*TEMP/(1.0D0+4.0D0*TEMP**2)**16 |
---|
| 1094 | GO TO 330 |
---|
| 1095 | 310 YP(1) = -2.0D0*Y(1) |
---|
| 1096 | GO TO 330 |
---|
| 1097 | ! 55: |
---|
| 1098 | ! PROBLEM F5 |
---|
| 1099 | 320 YP(1) = Y(1)*(4.0D0/(3.0D0*C2))*(SIGN(ABS(X- & |
---|
| 1100 | 1.0D0)**EX,X-1.0D0)+SIGN(ABS(X-2.0D0)**EX,X-2.0D0)+SIGN(ABS(X- & |
---|
| 1101 | 3.0D0)**EX,X-3.0D0)+SIGN(ABS(X-4.0D0)**EX,X-4.0D0)+SIGN(ABS(X- & |
---|
| 1102 | 5.0D0)**EX,X-5.0D0)+SIGN(ABS(X-6.0D0)**EX,X-6.0D0)+SIGN(ABS(X- & |
---|
| 1103 | 7.0D0)**EX,X-7.0D0)+SIGN(ABS(X-8.0D0)**EX,X-8.0D0)+SIGN(ABS(X- & |
---|
| 1104 | 9.0D0)**EX,X-9.0D0)+SIGN(ABS(X-10.0D0)**EX,X-10.0D0)+SIGN(ABS(X- & |
---|
| 1105 | 11.0D0)**EX,X-11.0D0)+SIGN(ABS(X-12.0D0)**EX,X-12.0D0)+SIGN(ABS(X- & |
---|
| 1106 | 13.0D0)**EX,X-13.0D0)+SIGN(ABS(X-14.0D0)**EX,X-14.0D0)+SIGN(ABS(X- & |
---|
| 1107 | 15.0D0)**EX,X-15.0D0)+SIGN(ABS(X-16.0D0)**EX,X-16.0D0)+SIGN(ABS(X- & |
---|
| 1108 | 17.0D0)**EX,X-17.0D0)+SIGN(ABS(X-18.0D0)**EX,X-18.0D0)+SIGN(ABS(X- & |
---|
| 1109 | 19.0D0)**EX,X-19.0D0)) |
---|
| 1110 | 330 CONTINUE |
---|
| 1111 | IF (IWT<0) GO TO 340 |
---|
| 1112 | DO I = 1, N |
---|
| 1113 | YP(I) = YP(I)/W(I) |
---|
| 1114 | Y(I) = YTEMP(I) |
---|
| 1115 | END DO |
---|
| 1116 | 340 CONTINUE |
---|
| 1117 | RETURN |
---|
| 1118 | END SUBROUTINE FCN |
---|
| 1119 | |
---|
| 1120 | END MODULE NONSTIFFSET |
---|
| 1121 | |
---|
| 1122 | !****************************************************************** |
---|
| 1123 | |
---|
| 1124 | PROGRAM DEMONONSTIFF |
---|
| 1125 | |
---|
| 1126 | USE NONSTIFFSET |
---|
| 1127 | USE DVODE_F90_M |
---|
| 1128 | |
---|
| 1129 | IMPLICIT NONE |
---|
| 1130 | INTEGER ITASK, ISTATE, ISTATS, NEQ, I, CLASS, PROBLEM, MYID, ITEST, ITOL |
---|
| 1131 | DOUBLE PRECISION RSTATS, T, TOUT, HBEGIN, HBOUND, TBEGIN, TEND, Y, EPS, & |
---|
| 1132 | YINIT, YFINAL, RELERR_TOLERANCES, ABSERR_TOLERANCES, AERROR |
---|
| 1133 | LOGICAL USEW, USEHBEGIN, ERRORS |
---|
| 1134 | DIMENSION Y(51), RSTATS(22), ISTATS(31), YINIT(51), YFINAL(51), MYID(55) |
---|
| 1135 | DIMENSION RELERR_TOLERANCES(51), ABSERR_TOLERANCES(51), AERROR(51) |
---|
| 1136 | TYPE (VODE_OPTS) :: OPTIONS |
---|
| 1137 | DATA MYID/1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 11, 12, 13, 14, 15, 0, 0, 0, 0, & |
---|
| 1138 | 0, 21, 22, 23, 24, 25, 0, 0, 0, 0, 0, 31, 32, 33, 34, 35, 0, 0, 0, 0, & |
---|
| 1139 | 0, 41, 42, 43, 44, 45, 0, 0, 0, 0, 0, 51, 52, 53, 54, 55/ |
---|
| 1140 | |
---|
| 1141 | OPEN (UNIT=6,FILE='nonstiffoptions.dat') |
---|
| 1142 | |
---|
| 1143 | ERRORS = .FALSE. |
---|
| 1144 | |
---|
| 1145 | DO ITOL = 4, 10, 2 |
---|
| 1146 | WRITE(6,*) 'Results for ITOL = ', ITOL, ' follow.' |
---|
| 1147 | DO ITEST = 1, 55 |
---|
| 1148 | ID = MYID(ITEST) |
---|
| 1149 | IF (ID==0) GO TO 20 |
---|
| 1150 | WRITE (6,90010) |
---|
| 1151 | CLASS = ID/10 |
---|
| 1152 | PROBLEM = ID - 10*CLASS |
---|
| 1153 | IF (CLASS==0) THEN |
---|
| 1154 | WRITE (6,90000) PROBLEM |
---|
| 1155 | ELSE IF (CLASS==1) THEN |
---|
| 1156 | WRITE (6,90001) PROBLEM |
---|
| 1157 | ELSE IF (CLASS==2) THEN |
---|
| 1158 | WRITE (6,90002) PROBLEM |
---|
| 1159 | ELSE IF (CLASS==3) THEN |
---|
| 1160 | WRITE (6,90003) PROBLEM |
---|
| 1161 | ELSE IF (CLASS==4) THEN |
---|
| 1162 | WRITE (6,90004) PROBLEM |
---|
| 1163 | ELSE IF (CLASS==5) THEN |
---|
| 1164 | WRITE (6,90005) PROBLEM |
---|
| 1165 | END IF |
---|
| 1166 | ! Scale the odes? |
---|
| 1167 | USEW = .TRUE. |
---|
| 1168 | ! Use the IVALU starting step size? |
---|
| 1169 | USEHBEGIN = .TRUE. |
---|
| 1170 | IWT = -1 |
---|
| 1171 | IF (USEW) IWT = 1 |
---|
| 1172 | CALL IVALU(TBEGIN,TEND,HBEGIN,HBOUND,YINIT) |
---|
| 1173 | IF (.NOT. USEHBEGIN) HBEGIN = 0.0D0 |
---|
| 1174 | NEQ = N |
---|
| 1175 | T = TBEGIN |
---|
| 1176 | TOUT = TEND |
---|
| 1177 | Y(1:NEQ) = YINIT(1:NEQ) |
---|
| 1178 | EPS = 1.0D0 / 10.0D0**ITOL |
---|
| 1179 | RELERR_TOLERANCES(1:NEQ) = EPS |
---|
| 1180 | ABSERR_TOLERANCES(1:NEQ) = EPS |
---|
| 1181 | WRITE (6,90007) ID, TBEGIN, TEND, HBEGIN, HBOUND, IWT, N, EPS, Y(1:NEQ) |
---|
| 1182 | ITASK = 1 |
---|
| 1183 | ISTATE = 1 |
---|
| 1184 | OPTIONS = SET_OPTS(RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ), & |
---|
| 1185 | ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),MXSTEP=100000,H0=HBEGIN, & |
---|
| 1186 | HMAX=HBOUND) |
---|
| 1187 | 10 CONTINUE |
---|
| 1188 | CALL DVODE_F90(DERIVS,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS) |
---|
| 1189 | IF (ISTATE<0) THEN |
---|
| 1190 | WRITE (6,90006) ISTATE |
---|
| 1191 | ERRORS = .TRUE. |
---|
| 1192 | ! STOP |
---|
| 1193 | END IF |
---|
| 1194 | CALL GET_STATS(RSTATS,ISTATS) |
---|
| 1195 | WRITE (6,90009) ISTATS(11), ISTATS(12) |
---|
| 1196 | IF (TOUT<TEND) GO TO 10 |
---|
| 1197 | CALL EVALU(YFINAL) |
---|
| 1198 | DO I = 1, NEQ |
---|
| 1199 | AERROR(I) = ABS(Y(I)-YFINAL(I)) |
---|
| 1200 | END DO |
---|
| 1201 | WRITE (6,90008) (I,Y(I),YFINAL(I),AERROR(I),I=1,NEQ) |
---|
| 1202 | 20 END DO ! End of ITEST Loop |
---|
| 1203 | END DO ! End of ITOL Loop |
---|
| 1204 | |
---|
| 1205 | IF (ERRORS) THEN |
---|
| 1206 | WRITE(6,90011) |
---|
| 1207 | ELSE |
---|
| 1208 | WRITE(6,90012) |
---|
| 1209 | END IF |
---|
| 1210 | |
---|
| 1211 | ! Format statements for this problem: |
---|
| 1212 | 90000 FORMAT (' Class/Problem = A',I1) |
---|
| 1213 | 90001 FORMAT (' Class/Problem = B',I1) |
---|
| 1214 | 90002 FORMAT (' Class/Problem = C',I1) |
---|
| 1215 | 90003 FORMAT (' Class/Problem = D',I1) |
---|
| 1216 | 90004 FORMAT (' Class/Problem = E',I1) |
---|
| 1217 | 90005 FORMAT (' Class/Problem = F',I1) |
---|
| 1218 | 90006 FORMAT (' An error occurred in VODE_F90. ISTATE = ',I3) |
---|
| 1219 | 90007 FORMAT (' Problem ID = ',I3,/,' Initial time = ',D15.5,/, & |
---|
| 1220 | ' Final time = ',D15.5,/,' Initial stepsize = ',D15.5,/, & |
---|
| 1221 | ' Maximum stepsize = ',D15.5,/,' IWT flag = ',I3,/, & |
---|
| 1222 | ' Number of odes = ',I3,/,' Error tolerance = ',D15.5,/, & |
---|
| 1223 | ' Initial solution = ',/,(D15.5)) |
---|
| 1224 | 90008 FORMAT (' Computed and reference solutions and absolute', & |
---|
| 1225 | ' errors follow:',/,(I3,3D15.5)) |
---|
| 1226 | 90009 FORMAT (' Steps = ',I10,' f-s = ',I10) |
---|
| 1227 | 90010 FORMAT (' _________________________________________') |
---|
| 1228 | 90011 FORMAT (' Errors occurred.') |
---|
| 1229 | 90012 FORMAT (' No errors occurred.') |
---|
| 1230 | STOP |
---|
| 1231 | |
---|
| 1232 | END PROGRAM DEMONONSTIFF |
---|