mmid.f

55 lines · 1 definitions · on GitHub (7be4133)

Python version: mmid.py, 72 lines, on the right (switch to Python above the code).

Purpose

The modified midpoint method of Numerical Recipes ("this SR from Numerical Recipes, p.562"): it advances three variables $y$ from $x$ to $x + H$ in $n$ substeps. It is the inner step of the Bulirsch–Stoer integrator BSSTEP, which ENVINT uses to integrate the envelope.

How it fits

Who calls it: only BSSTEP, with a sequence of $n$ (NSEQ).
What it calls: the derivative routine passed in as DERIV (an EXTERNAL argument; the site cannot resolve it statically). All the arguments after DERIV (B, FPL, FTL, GL, LATMO, …, KSAHA) are only passed through to it.
Output: YOUT, the values at $x+H$.

The method

With $h = H/n$, $z_0 = y(x)$, $z_1 = z_0 + h\,y'(x)$, $$z_{m+1} = z_{m-1} + 2h\,f(x + m h, z_m),\quad m = 1,\dots,n-1,\qquad y(x+H) \approx \tfrac12\left[z_n + z_{n-1} + h\,f(x+H, z_n)\right].$$ YOUT holds the derivatives in between, to save an array.

Design

The size is fixed at three variables (Y(3)), whatever NVAR is.

Definitions

subroutineMMID

lines 8–55 · GitHub

Called by BSSTEP

External calls (1)

DERIV ×2

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE.

1 2 3C 4C 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6C MMID 7C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 8 SUBROUTINE MMID(Y,DYDX,NVAR,XS,HTOT,NSTEP,YOUT,DERIV,B,FPL,FTL, 9 * * GL,LATMO,LDERIV,LOCOND,LSAVE,RL,TEFFL,X,Z,KOUNT,KSAHA) 10C 10/20/26 MHP via Claude 11C code changed to replace implicit typing with IMPLICIT NONE and 12C explicit declarations. 13 IMPLICIT NONE 14 REAL*8 B,DYDX,FPL,FTL,GL,H,H2,HTOT,RL,SWAP,TEFFL,X,X0,XS,Y,YM,YN, 15 * * YOUT,Z 16 INTEGER I,KOUNT,KSAHA,N,NSTEP,NVAR 17 LOGICAL*4 LATMO,LDERIV,LOCOND,LSAVE 18 19 DIMENSION Y(3),DYDX(3),YOUT(3),YM(3),YN(3) 20 EXTERNAL DERIV 21 SAVE 22 23C INTEGRATES THE DEPENDENT VARIABLES Y FROM XS TO XS + HTOT IN 24C NSTEP INCREMENTS. INPUT ARE THE Y'S AND DY/DX'S AT XS; VALUES OF 25C Y AT XS + HTOT ARE STORED IN YOUT. DERIVATIVES ARE CALCULATED IN 26C SUBROUTINE DERIV. THIS SR FROM NUMERICAL RECIPES, P.562. 27C H IS THE SIZE OF EACH SMALL STEP. 28 H = HTOT/DFLOAT(NSTEP) 29C FIRST STEP 30 DO 10 I = 1,NVAR 31 YM(I) = Y(I) 32 YN(I) = Y(I) + DYDX(I)*H 33 10 10 CONTINUE 34 X0 = XS + H 35C YOUT TEMPORARILY USED FOR STORAGE OF DERIVATIVES. 36 CALL DERIV(X0,YN,YOUT,B,FPL,FTL,GL,LATMO,LDERIV,LOCOND,LSAVE, 37 * * RL,TEFFL,X,Z,KOUNT,KSAHA) 38 H2 = 2.0D0*H 39C GENERAL STEP. 40 DO 30 N = 2,NSTEP 41 DO 20 I = 1,NVAR 42 SWAP = YM(I) + H2*YOUT(I) 43 YM(I) = YN(I) 44 YN(I) = SWAP 45 20 20 CONTINUE 46 X0 = X0 + H 47 CALL DERIV(X0,YN,YOUT,B,FPL,FTL,GL,LATMO,LDERIV,LOCOND,LSAVE, 48 * * RL,TEFFL,X,Z,KOUNT,KSAHA) 49 30 30 CONTINUE 50C LAST STEP. 51 DO 40 I = 1,NVAR 52 YOUT(I) = 0.5D0*(YM(I) + YN(I) + H*YOUT(I)) 53 40 40 CONTINUE 54 RETURN 55 END

1"""Python version of mmid.f (YREC 7be4133): the modified midpoint step of 2Numerical Recipes, the inner step of BSSTEP.""" 3 4from collections.abc import Callable 5 6import numpy as np 7 8# The Fortran's local arrays, kept by SAVE (Fortran line 19). 9ym = np.zeros(4) 10yn = np.zeros(4) 11 12 13def _call_deriv(deriv, x0, y, dydx, b, fpl, ftl, gl, latmo, lderiv, locond, 14 lsave, rl, teffl, x, z, kount, ksaha): 15 # DERIV is QENV, which returns (rl, kenv, ksaha), or QATM, which returns 16 # (katm, ksaha); the Fortran passes RL, KOUNT and KSAHA by reference. 17 result = deriv(x0, y, dydx, b, fpl, ftl, gl, latmo, lderiv, locond, 18 lsave, rl, teffl, x, z, kount, ksaha) 19 if len(result) == 3: 20 rl, kount, ksaha = result 21 else: 22 kount, ksaha = result 23 return rl, kount, ksaha 24 25 26def mmid(y: np.ndarray, dydx: np.ndarray, nvar: int, xs: float, htot: float, nstep: int, yout: np.ndarray, 27 deriv: Callable, b: float, fpl: float, ftl: float, gl: float, latmo: bool, lderiv: bool, locond: bool, 28 lsave: bool, rl: float, teffl: float, x: float, z: float, kount: int, ksaha: int) -> tuple[float, int, int]: 29 """Advance y from xs to xs + htot in nstep substeps by the modified 30 midpoint method; the result goes to yout. 31 32 Translates: mmid.f#MMID (7be4133) 33 34 Arguments: y and dydx at xs (three variables, index 0 unused; nvar of 35 them used); the step htot and the number of substeps nstep; yout, set 36 in place; deriv, the derivative routine (qenv or qatm); the arguments 37 from b to ksaha are only passed on to deriv. 38 39 Returns (rl, kount, ksaha), as deriv leaves them. 40 41 """ 42 43 # The substep (Fortran line 28). 44 h = htot / float(nstep) 45 46 # First step (Fortran lines 29–34). 47 for i in range(1, nvar + 1): # loop 10 48 ym[i] = y[i] 49 yn[i] = y[i] + dydx[i] * h 50 x0 = xs + h 51 52 # yout holds the derivatives in between (Fortran lines 35–38). 53 rl, kount, ksaha = _call_deriv(deriv, x0, yn, yout, b, fpl, ftl, gl, 54 latmo, lderiv, locond, lsave, rl, teffl, 55 x, z, kount, ksaha) 56 h2 = 2.0 * h 57 58 # General step (Fortran lines 39–49). 59 for n in range(2, nstep + 1): # loop 30 60 for i in range(1, nvar + 1): 61 swap = ym[i] + h2 * yout[i] 62 ym[i] = yn[i] 63 yn[i] = swap 64 x0 = x0 + h 65 rl, kount, ksaha = _call_deriv(deriv, x0, yn, yout, b, fpl, ftl, gl, 66 latmo, lderiv, locond, lsave, rl, 67 teffl, x, z, kount, ksaha) 68 69 # Last step (Fortran lines 50–54). 70 for i in range(1, nvar + 1): # loop 40 71 yout[i] = 0.5 * (ym[i] + yn[i] + h * yout[i]) 72 return rl, kount, ksaha