1*DECK SINTRP 2 SUBROUTINE SINTRP (X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, IVC, 3 + IV, KGI, GI, ALPHA, OG, OW, OX, OY) 4C***BEGIN PROLOGUE SINTRP 5C***PURPOSE Approximate the solution at XOUT by evaluating the 6C polynomial computed in STEPS at XOUT. Must be used in 7C conjunction with STEPS. 8C***LIBRARY SLATEC (DEPAC) 9C***CATEGORY I1A1B 10C***TYPE SINGLE PRECISION (SINTRP-S, DINTP-D) 11C***KEYWORDS ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE, 12C ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR, 13C SMOOTH INTERPOLANT 14C***AUTHOR Watts, H. A., (SNLA) 15C***DESCRIPTION 16C 17C The methods in subroutine STEPS approximate the solution near X 18C by a polynomial. Subroutine SINTRP approximates the solution at 19C XOUT by evaluating the polynomial there. Information defining this 20C polynomial is passed from STEPS so SINTRP cannot be used alone. 21C 22C Subroutine STEPS is completely explained and documented in the text, 23C "Computer Solution of Ordinary Differential Equations, the Initial 24C Value Problem" by L. F. Shampine and M. K. Gordon. 25C 26C Input to SINTRP -- 27C 28C The user provides storage in the calling program for the arrays in 29C the call list 30C DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN) 31C AND ALPHA(12),OG(13),OW(12),GI(11),IV(10) 32C and defines 33C XOUT -- point at which solution is desired. 34C The remaining parameters are defined in STEPS and passed to 35C SINTRP from that subroutine 36C 37C Output from SINTRP -- 38C 39C YOUT(*) -- solution at XOUT 40C YPOUT(*) -- derivative of solution at XOUT 41C The remaining parameters are returned unaltered from their input 42C values. Integration with STEPS may be continued. 43C 44C***REFERENCES H. A. Watts, A smoother interpolant for DE/STEP, INTRP 45C II, Report SAND84-0293, Sandia Laboratories, 1984. 46C***ROUTINES CALLED (NONE) 47C***REVISION HISTORY (YYMMDD) 48C 840201 DATE WRITTEN 49C 890831 Modified array declarations. (WRB) 50C 890831 REVISION DATE from Version 3.2 51C 891214 Prologue converted to Version 4.0 format. (BAB) 52C 920501 Reformatted the REFERENCES section. (WRB) 53C***END PROLOGUE SINTRP 54C 55 DIMENSION Y(*),YOUT(*),YPOUT(*),PHI(NEQN,16),OY(*) 56 DIMENSION G(13),C(13),W(13),OG(13),OW(12),ALPHA(12),GI(11),IV(10) 57C 58C***FIRST EXECUTABLE STATEMENT SINTRP 59 KP1 = KOLD + 1 60 KP2 = KOLD + 2 61C 62 HI = XOUT - OX 63 H = X - OX 64 XI = HI/H 65 XIM1 = XI - 1. 66C 67C INITIALIZE W(*) FOR COMPUTING G(*) 68C 69 XIQ = XI 70 DO 10 IQ = 1,KP1 71 XIQ = XI*XIQ 72 TEMP1 = IQ*(IQ+1) 73 10 W(IQ) = XIQ/TEMP1 74C 75C COMPUTE THE DOUBLE INTEGRAL TERM GDI 76C 77 IF (KOLD .LE. KGI) GO TO 50 78 IF (IVC .GT. 0) GO TO 20 79 GDI = 1.0/TEMP1 80 M = 2 81 GO TO 30 82 20 IW = IV(IVC) 83 GDI = OW(IW) 84 M = KOLD - IW + 3 85 30 IF (M .GT. KOLD) GO TO 60 86 DO 40 I = M,KOLD 87 40 GDI = OW(KP2-I) - ALPHA(I)*GDI 88 GO TO 60 89 50 GDI = GI(KOLD) 90C 91C COMPUTE G(*) AND C(*) 92C 93 60 G(1) = XI 94 G(2) = 0.5*XI*XI 95 C(1) = 1.0 96 C(2) = XI 97 IF (KOLD .LT. 2) GO TO 90 98 DO 80 I = 2,KOLD 99 ALP = ALPHA(I) 100 GAMMA = 1.0 + XIM1*ALP 101 L = KP2 - I 102 DO 70 JQ = 1,L 103 70 W(JQ) = GAMMA*W(JQ) - ALP*W(JQ+1) 104 G(I+1) = W(1) 105 80 C(I+1) = GAMMA*C(I) 106C 107C DEFINE INTERPOLATION PARAMETERS 108C 109 90 SIGMA = (W(2) - XIM1*W(1))/GDI 110 RMU = XIM1*C(KP1)/GDI 111 HMU = RMU/H 112C 113C INTERPOLATE FOR THE SOLUTION -- YOUT 114C AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT 115C 116 DO 100 L = 1,NEQN 117 YOUT(L) = 0.0 118 100 YPOUT(L) = 0.0 119 DO 120 J = 1,KOLD 120 I = KP2 - J 121 GDIF = OG(I) - OG(I-1) 122 TEMP2 = (G(I) - G(I-1)) - SIGMA*GDIF 123 TEMP3 = (C(I) - C(I-1)) + RMU*GDIF 124 DO 110 L = 1,NEQN 125 YOUT(L) = YOUT(L) + TEMP2*PHI(L,I) 126 110 YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I) 127 120 CONTINUE 128 DO 130 L = 1,NEQN 129 YOUT(L) = ((1.0 - SIGMA)*OY(L) + SIGMA*Y(L)) + 130 1 H*(YOUT(L) + (G(1) - SIGMA*OG(1))*PHI(L,1)) 131 130 YPOUT(L) = HMU*(OY(L) - Y(L)) + 132 1 (YPOUT(L) + (C(1) + RMU*OG(1))*PHI(L,1)) 133C 134 RETURN 135 END 136