home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / jpl_c / rknstep.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-05-18  |  2.1 KB  |  62 lines
  1. /* 1.0  11-19-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************
  7.  *    Integration step for N first-order ordinary differential equations
  8.  *    using the 4th-order Runge-Kutta formulas, given in
  9.  *
  10.  *    HANDBOOK OF MATHEMATICAL FUNCTIONS, U. S. Dept. of Commerce,
  11.  *    National Bureau of Standards Applied Mathematics Series No. 55,
  12.  *    formula 25.5.10, page 896, June, 1964,
  13.  *
  14.  *    on each equation y'[i] = f(i, t, y[i]), i = 0, ... , N-1, as time
  15.  *    progresses.
  16.  *
  17.  *----------------------------------------------------------------------*/
  18.  
  19. #include "defs.h"
  20. #include "stdtyp.h"
  21.  
  22. /************************************************************************/
  23.     VOID
  24. rknstep(N, f, h, t, y1, y2, dy, k, ys)
  25.         /* Store y2[i] = y1[i](t+h), for y1[i]' = f(i, t, y1),
  26.            using 4th-order Runge-Kutta method, i = 0, ... , N.    */
  27. /*----------------------------------------------------------------------*/
  28. double (*f)(), h, t, y1[], y2[], dy[], k[], ys[];
  29. /*
  30.  *    Note: dy[] is a scratch array needed to hold intermediate 
  31.  *    calculations accumulating the dy value; k[] holds coefficients
  32.  *    k1, k2, k3, and k4, successively; and the ys[] is for y approximants,
  33.  *    for each i.
  34.  */
  35. {
  36.     int i;
  37.     double h2, ldexp();
  38.  
  39.     h2 = ldexp(h, -1);
  40.     for (i = 0; i < N; i++)        /* i is equation index.    */
  41.     {    dy[i] = h * (*f)(i, t, y1);        /* k1    */
  42.         ys[i] = y1[i] + dy[i] / 2;
  43.     }
  44.     for (i = 0; i < N; i++)
  45.     {    k[i] = h * (*f)(i, t + h2, ys);        /* k2    */
  46.         dy[i] += 2 * k[i];
  47.     }
  48.     for (i = 0; i < N; i++)
  49.         ys[i] = y1[i] + k[i] / 2.0;
  50.     for (i = 0; i < N; i++)
  51.     {    k[i] = h * (*f)(i, t + h2, ys);        /* k3    */
  52.         dy[i] += 2 * k[i];
  53.     }
  54.     for (i = 0; i < N; i++)
  55.         ys[i] = y1[i] + k[i];
  56.     for (i = 0; i < N; i++)
  57.         dy[i] += h * (*f)(i, t + h, ys);    /* k4    */
  58.                 /* dy is (k1 + 2 k2 + 2 k3 + k4) */
  59.     for (i = 0; i < N; i++)
  60.         y2[i] = y1[i] + dy[i] / 6.0;
  61. }
  62.