GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
sintdy.f
Go to the documentation of this file.
1  SUBROUTINE sintdy (T, K, YH, NYH, DKY, IFLAG)
2 C***BEGIN PROLOGUE SINTDY
3 C***SUBSIDIARY
4 C***PURPOSE Interpolate solution derivatives.
5 C***TYPE SINGLE PRECISION (SINTDY-S, DINTDY-D)
6 C***AUTHOR Hindmarsh, Alan C., (LLNL)
7 C***DESCRIPTION
8 C
9 C SINTDY computes interpolated values of the K-th derivative of the
10 C dependent variable vector y, and stores it in DKY. This routine
11 C is called within the package with K = 0 and T = TOUT, but may
12 C also be called by the user for any K up to the current order.
13 C (See detailed instructions in the usage documentation.)
14 C
15 C The computed values in DKY are gotten by interpolation using the
16 C Nordsieck history array YH. This array corresponds uniquely to a
17 C vector-valued polynomial of degree NQCUR or less, and DKY is set
18 C to the K-th derivative of this polynomial at T.
19 C The formula for DKY is:
20 C q
21 C DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
22 C j=K
23 C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
24 C The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
25 C communicated by COMMON. The above sum is done in reverse order.
26 C IFLAG is returned negative if either K or T is out of bounds.
27 C
28 C***SEE ALSO SLSODE
29 C***ROUTINES CALLED XERRWV
30 C***COMMON BLOCKS SLS001
31 C***REVISION HISTORY (YYMMDD)
32 C 791129 DATE WRITTEN
33 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
34 C 890503 Minor cosmetic changes. (FNF)
35 C 930809 Renamed to allow single/double precision versions. (ACH)
36 C 010412 Reduced size of Common block /SLS001/. (ACH)
37 C 031105 Restored 'own' variables to Common block /SLS001/, to
38 C enable interrupt/restart feature. (ACH)
39 C 050427 Corrected roundoff decrement in TP. (ACH)
40 C***END PROLOGUE SINTDY
41 C**End
42  INTEGER k, nyh, iflag
43  REAL t, yh, dky
44  dimension yh(nyh,*), dky(*)
45  INTEGER iownd, iowns,
46  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
47  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
48  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
49  REAL rowns,
50  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
51  COMMON /sls001/ rowns(209),
52  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
53  2 iownd(6), iowns(6),
54  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
55  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
56  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
57  INTEGER i, ic, j, jb, jb2, jj, jj1, jp1
58  REAL c, r, s, tp
59  CHARACTER*80 msg
60 C
61 C***FIRST EXECUTABLE STATEMENT SINTDY
62  iflag = 0
63  IF (k .LT. 0 .OR. k .GT. nq) go to 80
64  tp = tn - hu - 100.0e0*uround*sign(abs(tn) + abs(hu), hu)
65  IF ((t-tp)*(t-tn) .GT. 0.0e0) go to 90
66 C
67  s = (t - tn)/h
68  ic = 1
69  IF (k .EQ. 0) go to 15
70  jj1 = l - k
71  DO 10 jj = jj1,nq
72  10 ic = ic*jj
73  15 c = ic
74  DO 20 i = 1,n
75  20 dky(i) = c*yh(i,l)
76  IF (k .EQ. nq) go to 55
77  jb2 = nq - k
78  DO 50 jb = 1,jb2
79  j = nq - jb
80  jp1 = j + 1
81  ic = 1
82  IF (k .EQ. 0) go to 35
83  jj1 = jp1 - k
84  DO 30 jj = jj1,j
85  30 ic = ic*jj
86  35 c = ic
87  DO 40 i = 1,n
88  40 dky(i) = c*yh(i,jp1) + s*dky(i)
89  50 CONTINUE
90  IF (k .EQ. 0) RETURN
91  55 r = h**(-k)
92  DO 60 i = 1,n
93  60 dky(i) = r*dky(i)
94  RETURN
95 C
96  80 CALL xerrwd('SINTDY- K (=I1) illegal ',
97  1 30, 51, 0, 1, k, 0, 0, 0.0e0, 0.0e0)
98  iflag = -1
99  RETURN
100  90 CALL xerrwd('SINTDY- T (=R1) illegal ',
101  1 30, 52, 0, 0, 0, 0, 1, t, 0.0e0)
102  CALL xerrwd(
103  1 ' T not in interval TCUR - HU (= R1) to TCUR (=R2) ',
104  1 60, 52, 0, 0, 0, 0, 2, tp, tn)
105  iflag = -2
106  RETURN
107 C----------------------- END OF SUBROUTINE SINTDY ----------------------
108  END