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
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 iownd, iowns,
5  1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
6  2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
7  INTEGER i, ic, j, jb, jb2, jj, jj1, jp1
8  DOUBLE PRECISION t, yh, dky
9  DOUBLE PRECISION rowns,
10  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
11  DOUBLE PRECISION c, r, s, tp
12  dimension yh(nyh,*), dky(*)
13  COMMON /ls0001/ rowns(209),
14  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
15  3 iownd(14), iowns(6),
16  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
17  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
18 C-----------------------------------------------------------------------
19 C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE
20 C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY. THIS ROUTINE
21 C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY
22 C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER.
23 C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.)
24 C-----------------------------------------------------------------------
25 C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE
26 C NORDSIECK HISTORY ARRAY YH. THIS ARRAY CORRESPONDS UNIQUELY TO A
27 C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET
28 C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T.
29 C THE FORMULA FOR DKY IS..
30 C Q
31 C DKY(I) = SUM C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
32 C J=K
33 C WHERE C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
34 C THE QUANTITIES NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE
35 C COMMUNICATED BY COMMON. THE ABOVE SUM IS DONE IN REVERSE ORDER.
36 C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS.
37 C-----------------------------------------------------------------------
38  iflag = 0
39  IF (k .LT. 0 .OR. k .GT. nq) go to 80
40  tp = tn - hu - 100.0d0*uround*(tn + hu)
41  IF ((t-tp)*(t-tn) .GT. 0.0d0) go to 90
42 C
43  s = (t - tn)/h
44  ic = 1
45  IF (k .EQ. 0) go to 15
46  jj1 = l - k
47  DO 10 jj = jj1,nq
48  10 ic = ic*jj
49  15 c = dble(ic)
50  DO 20 i = 1,n
51  20 dky(i) = c*yh(i,l)
52  IF (k .EQ. nq) go to 55
53  jb2 = nq - k
54  DO 50 jb = 1,jb2
55  j = nq - jb
56  jp1 = j + 1
57  ic = 1
58  IF (k .EQ. 0) go to 35
59  jj1 = jp1 - k
60  DO 30 jj = jj1,j
61  30 ic = ic*jj
62  35 c = dble(ic)
63  DO 40 i = 1,n
64  40 dky(i) = c*yh(i,jp1) + s*dky(i)
65  50 CONTINUE
66  IF (k .EQ. 0) RETURN
67  55 r = h**(-k)
68  DO 60 i = 1,n
69  60 dky(i) = r*dky(i)
70  RETURN
71 C
72  80 CALL xerrwd('INTDY-- K (=I1) ILLEGAL ',
73  1 30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
74  iflag = -1
75  RETURN
76  90 CALL xerrwd('INTDY-- T (=R1) ILLEGAL ',
77  1 30, 52, 0, 0, 0, 0, 1, t, 0.0d0)
78  CALL xerrwd(
79  1 ' T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2) ',
80  1 60, 52, 0, 0, 0, 0, 2, tp, tn)
81  iflag = -2
82  RETURN
83 C----------------------- END OF SUBROUTINE INTDY -----------------------
84  END