1 SUBROUTINE INTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,PSI) 2C***BEGIN PROLOGUE INTRP 3C***DATE WRITTEN 740101 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***CATEGORY NO. Z 6C***KEYWORDS ADAMS METHOD,DEPAC,INITIAL VALUE PROBLEMS,ODE, 7C ORDINARY DIFFERENTIAL EQUATIONS,PREDICTOR-CORRECTOR 8C***AUTHOR SHAMPINE, L. F., (SNLA) 9C GORDON, M. K., (SNLA) 10C***PURPOSE Approximates the solution at XOUT by evaluating the 11C polynomial computed in STEP2 at XOUT. Must be used in 12C conjunction with STEP2. 13C***DESCRIPTION 14C 15C Written by L. F. Shampine and M. K. Gordon 16C 17C Abstract 18C 19C 20C The methods in subroutine STEP2 approximate the solution near X 21C by a polynomial. Subroutine INTRP approximates the solution at 22C XOUT by evaluating the polynomial there. Information defining this 23C polynomial is passed from STEP2 so INTRP cannot be used alone. 24C 25C This code is completely explained and documented in the text, 26C "Computer Solution of Ordinary Differential Equations, the Initial 27C Value Problem" by L. F. Shampine and M. K. Gordon. 28C Further details on use of this code are available in "Solving 29C Ordinary Differential Equations with ODE, STEP, and INTRP", 30C by L. F. Shampine and M. K. Gordon, SLA-73-1060. 31C 32C Input to INTRP -- 33C 34C The user provides storage in the calling program for the arrays in 35C the call list 36C DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12) 37C and defines 38C XOUT -- point at which solution is desired. 39C The remaining parameters are defined in STEP2 and passed to 40C INTRP from that subroutine 41C 42C Output from INTRP -- 43C 44C YOUT(*) -- solution at XOUT 45C YPOUT(*) -- derivative of solution at XOUT 46C The remaining parameters are returned unaltered from their input 47C values. Integration with STEP2 may be continued. 48C***REFERENCES SHAMPINE, L. F., GORDON, M. K., *SOLVING 49C ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, 50C AND INTRP*, SLA-73-1060, SANDIA LABORATORIES, 51C 1973. 52C***ROUTINES CALLED (NONE) 53C***END PROLOGUE INTRP 54C 55C 56 DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12) 57 DIMENSION G(13),W(13),RHO(13) 58 DATA G(1)/1.0/,RHO(1)/1.0/ 59C 60C***FIRST EXECUTABLE STATEMENT INTRP 61 HI = XOUT - X 62 KI = KOLD + 1 63 KIP1 = KI + 1 64C 65C INITIALIZE W(*) FOR COMPUTING G(*) 66C 67 DO 5 I = 1,KI 68 TEMP1 = I 69 5 W(I) = 1.0/TEMP1 70 TERM = 0.0 71C 72C COMPUTE G(*) 73C 74 DO 15 J = 2,KI 75 JM1 = J - 1 76 PSIJM1 = PSI(JM1) 77 GAMMA = (HI + TERM)/PSIJM1 78 ETA = HI/PSIJM1 79 LIMIT1 = KIP1 - J 80 DO 10 I = 1,LIMIT1 81 10 W(I) = GAMMA*W(I) - ETA*W(I+1) 82 G(J) = W(1) 83 RHO(J) = GAMMA*RHO(JM1) 84 15 TERM = PSIJM1 85C 86C INTERPOLATE 87C 88 DO 20 L = 1,NEQN 89 YPOUT(L) = 0.0 90 20 YOUT(L) = 0.0 91 DO 30 J = 1,KI 92 I = KIP1 - J 93 TEMP2 = G(I) 94 TEMP3 = RHO(I) 95 DO 25 L = 1,NEQN 96 YOUT(L) = YOUT(L) + TEMP2*PHI(L,I) 97 25 YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I) 98 30 CONTINUE 99 DO 35 L = 1,NEQN 100 35 YOUT(L) = Y(L) + HI*YOUT(L) 101 RETURN 102 END 103