GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
intdy.f
Go to the documentation of this file.
1  SUBROUTINE intdy (T, K, YH, NYH, DKY, IFLAG)
2 CLLL. OPTIMIZE
3  INTEGER K, NYH, IFLAG
4  INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
5  1 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, CNYH,
6  2 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP
7  INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
8  2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
9  INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
10  DOUBLE PRECISION T, YH, DKY
11  DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
12  1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
13  DOUBLE PRECISION C, R, S, TP
14  dimension yh(nyh,*), dky(*)
15  COMMON /ls0001/ conit, crate, el(13), elco(13,12),
16  1 hold, rmax, tesco(3,12),
17  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
18  2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
19  3 mxstep, mxhnil, nhnil, ntrep, nslast, cnyh,
20  3 ialth, ipup, lmax, meo, nqnyh, nslp,
21  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
22  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
23 C-----------------------------------------------------------------------
24 C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE
25 C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY. THIS ROUTINE
26 C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY
27 C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER.
28 C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.)
29 C-----------------------------------------------------------------------
30 C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE
31 C NORDSIECK HISTORY ARRAY YH. THIS ARRAY CORRESPONDS UNIQUELY TO A
32 C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET
33 C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T.
34 C THE FORMULA FOR DKY IS..
35 C Q
36 C DKY(I) = SUM C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
37 C J=K
38 C WHERE C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
39 C THE QUANTITIES NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE
40 C COMMUNICATED BY COMMON. THE ABOVE SUM IS DONE IN REVERSE ORDER.
41 C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS.
42 C-----------------------------------------------------------------------
43  iflag = 0
44  IF (k .LT. 0 .OR. k .GT. nq) GO TO 80
45  tp = tn - hu - 100.0d0*uround*(tn + hu)
46  IF ((t-tp)*(t-tn) .GT. 0.0d0) GO TO 90
47 C
48  s = (t - tn)/h
49  ic = 1
50  IF (k .EQ. 0) GO TO 15
51  jj1 = l - k
52  DO 10 jj = jj1,nq
53  10 ic = ic*jj
54  15 c = dble(ic)
55  DO 20 i = 1,n
56  20 dky(i) = c*yh(i,l)
57  IF (k .EQ. nq) GO TO 55
58  jb2 = nq - k
59  DO 50 jb = 1,jb2
60  j = nq - jb
61  jp1 = j + 1
62  ic = 1
63  IF (k .EQ. 0) GO TO 35
64  jj1 = jp1 - k
65  DO 30 jj = jj1,j
66  30 ic = ic*jj
67  35 c = dble(ic)
68  DO 40 i = 1,n
69  40 dky(i) = c*yh(i,jp1) + s*dky(i)
70  50 CONTINUE
71  IF (k .EQ. 0) RETURN
72  55 r = h**(-k)
73  DO 60 i = 1,n
74  60 dky(i) = r*dky(i)
75  RETURN
76 C
77  80 CALL xerrwd('INTDY-- K (=I1) ILLEGAL ',
78  1 30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
79  iflag = -1
80  RETURN
81  90 CALL xerrwd('INTDY-- T (=R1) ILLEGAL ',
82  1 30, 52, 0, 0, 0, 0, 1, t, 0.0d0)
83  CALL xerrwd(
84  1 ' T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2) ',
85  1 60, 52, 0, 0, 0, 0, 2, tp, tn)
86  iflag = -2
87  RETURN
88 C----------------------- END OF SUBROUTINE INTDY -----------------------
89  END
subroutine intdy(T, K, YH, NYH, DKY, IFLAG)
Definition: intdy.f:2
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:4