GNU Octave  4.2.1
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
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 xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:3
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine intdy(T, K, YH, NYH, DKY, IFLAG)
Definition: intdy.f:1