Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/DVODE/nonstiffoptions.f90 @ 4685

Last change on this file since 4685 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk
File size: 33.8 KB
Line 
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
9410      IOUT = 6
95        WRITE (IOUT,FMT=90000) ID
96        STOP
97
98!     PROBLEM CLASS A
99
100!     1:
10120      CONTINUE
102!     PROBLEM A1
103        N = 1
104        W(1) = 0.100D+01
105        Y(1) = 1.0D0
106        GO TO 330
107!     2:
10830      CONTINUE
109!     PROBLEM A2
110        N = 1
111        W(1) = 0.100D+01
112        Y(1) = 1.0D0
113        GO TO 330
114!     3:
11540      CONTINUE
116!     PROBLEM A3
117        N = 1
118        W(1) = 0.271D+01
119        Y(1) = 1.0D0
120        GO TO 330
121!     4:
12250      CONTINUE
123!     PROBLEM A4
124        N = 1
125        W(1) = 0.177D+02
126        Y(1) = 1.0D0
127        GO TO 330
128!     5:
12960      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:
13970      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:
14880      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:
15990      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:
170100     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:
181110     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:
195120     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:
214130     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:
233140     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:
252150     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:
312160     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
379170     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
388180     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
393190     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
398200     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
403210     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
408220     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:
418230     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:
428240     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:
437250     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:
446260     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:
455270     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:
467280     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:
477290     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:
485300     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:
496310     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:
504320     CONTINUE
505!     PROBLEM F5
506        N = 1
507        W(1) = 0.100D+01
508        Y(1) = 1.0D0
509        HMAX = 20.0D0
510330     CONTINUE
511        IF (IWT<0.) GO TO 340
512        DO I = 1, N
513          Y(I) = Y(I)/W(I)
514        END DO
515340     CONTINUE
516        RETURN
517
51890000   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
56710      Y(1) = 2.061153353012535D-09
568        GO TO 310
569!     2:
570!        PROBLEM A2
57120      Y(1) = 2.182178902359887D-01
572        GO TO 310
573!     3:
574!        PROBLEM A3
57530      Y(1) = 2.491650271850414D+00
576        GO TO 310
577!     4:
578!        PROBLEM A4
57940      Y(1) = 1.773016648131483D+01
580        GO TO 310
581!     5:
582!        PROBLEM A5
58350      Y(1) = -7.887826688964196D-01
584        GO TO 310
585
586!     PROBLEM CLASS B
587
588!     11:
589!        PROBLEM B1
59060      Y(1) = 6.761876008576667D-01
591        Y(2) = 1.860816099640036D-01
592        GO TO 310
593!     12:
594!        PROBLEM B2
59570      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
60180      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
60790      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
613100     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
622110     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
635120     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
648130     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
661140     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
715150     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
751160     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
758170     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
765180     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
772190     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
779200     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
789210     Y(1) = 1.456723600728308D-01
790        Y(2) = -9.883500195574063D-02
791        GO TO 310
792!     42:
793!        PROBLEM E2
794220     Y(1) = 2.008149762174948D+00
795        Y(2) = -4.250887527320057D-02
796        GO TO 310
797!     43:
798!        PROBLEM E3
799230     Y(1) = -1.004178858647128D-01
800        Y(2) = 2.411400132095954D-01
801        GO TO 310
802!     44:
803!        PROBLEM E4
804240     Y(1) = 3.395091444646555D+01
805        Y(2) = 2.767822659672869D-01
806        GO TO 310
807!     45:
808!        PROBLEM E5
809250     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
817260     Y(1) = -1.294460621213470D1
818        Y(2) = -2.208575158908672D-15
819        GO TO 310
820!     52:
821!        PROBLEM F2
822270     Y(1) = 70.03731057008607D0
823        GO TO 310
824!     53:
825!        PROBLEM F3
826280     Y(1) = -3.726957553088175D-1
827        Y(2) = -6.230137949234190D-1
828        GO TO 310
829!     54:
830!        PROBLEM F4
831290     Y(1) = 9.815017249707434D-11
832        GO TO 310
833!     55:
834!        PROBLEM F5
835300     Y(1) = 1.0D0
836310     CONTINUE
837        IF (IWT<0) GO TO 320
838        DO I = 1, N
839          Y(I) = Y(I)/W(I)
840        END DO
841320     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
88810      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
89920      YP(1) = -Y(1)
900        GO TO 330
901!     2:
902!        PROBLEM A2
90330      YP(1) = -.5D0*Y(1)*Y(1)*Y(1)
904        GO TO 330
905!     3:
906!        PROBLEM A3
90740      YP(1) = Y(1)*COS(X)
908        GO TO 330
909!     4:
910!        PROBLEM A4
91150      YP(1) = (1.0D0-Y(1)/20.0D0)*Y(1)/4.D0
912        GO TO 330
913!     5:
914!        PROBLEM A5
91560      YP(1) = (Y(1)-X)/(Y(1)+X)
916        GO TO 330
917
918!     PROBLEM CLASS B
919
920!     11:
921!        PROBLEM B1
92270      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
92880      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
93490      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
941100     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
948110     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
957120     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
965130     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
973140     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
981150     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
989160     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
1000170         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)
1003180         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
1027190     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
1039200     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
1044210     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
1049220     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
1054230     YP(1) = Y(2)
1055        YP(2) = .032D0 - .4D0*Y(2)*Y(2)
1056        GO TO 330
1057!     45:
1058!        PROBLEM E5
1059240     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
1067250     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
1073260     YP(2) = YP(2) + 1.0D0
1074        GO TO 330
1075!     52:
1076!        PROBLEM F2
1077270     ITEMP = INT(X)
1078        IF ((ITEMP/2)*2==ITEMP) GO TO 280
1079        YP(1) = 55.0D0 - .50D0*Y(1)
1080        GO TO 330
1081280     YP(1) = 55.0D0 - 1.50D0*Y(1)
1082        GO TO 330
1083!     53:
1084!        PROBLEM F3
1085290     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
1091300     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
1095310     YP(1) = -2.0D0*Y(1)
1096        GO TO 330
1097!     55:
1098!        PROBLEM F5
1099320     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))
1110330     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
1116340     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)
118710      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)
120220    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:
121290000 FORMAT (' Class/Problem = A',I1)
121390001 FORMAT (' Class/Problem = B',I1)
121490002 FORMAT (' Class/Problem = C',I1)
121590003 FORMAT (' Class/Problem = D',I1)
121690004 FORMAT (' Class/Problem = E',I1)
121790005 FORMAT (' Class/Problem = F',I1)
121890006 FORMAT (' An error occurred in VODE_F90. ISTATE = ',I3)
121990007 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))
122490008 FORMAT (' Computed and reference solutions and absolute', &
1225        ' errors follow:',/,(I3,3D15.5))
122690009 FORMAT (' Steps = ',I10,' f-s = ',I10)
122790010 FORMAT (' _________________________________________')
122890011 FORMAT (' Errors occurred.')
122990012 FORMAT (' No errors occurred.')
1230      STOP
1231
1232    END PROGRAM DEMONONSTIFF
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |