4.7.3 $B

($B$+$J$j$$$$2C8:$J%W%m%0%i%`$G$"$k$,!"(B $B$"$($F$3$3$K4^$a$kM}M3$N0l$D$O!"(B $B8x<0Cf$N78?t$K=q$-4V0c$$$,$J$$$3$H$r3NG'$G$-$k$h$&$K$9$k$?$a$G$"$k!#(B)

/*
 * rkf.c --- RKF45 $B$N%5%s%W%k(B ($B$^$@L$40@.$G$9(B)
 */

#include <math.h>

#define RKF45_STAGE 6

static double rkf45_alpha[RKF45_STAGE] = {
  0.0, 1.0/4, 3.0/8, 12.0/13, 1.0, 1.0/2};

static double rkf45_beta[RKF45_STAGE][RKF45_STAGE-1] = {
  {   0.0,          0.0,                      0.0,        0.0},
  {   1.0/4,        0.0,                      0.0,        0.0},
  {   3.0/32,       9.0/32,                   0.0,        0.0},
  {1932.0/2197, -7200.0/2197, 7296.0/2197,    0.0,        0.0},
  { 439.0/216,     -8.0,      3680.0/513,  -845.0/4104,   0.0},
  {  -8.0/27,       2.0,     -3544.0/2565, 1859.0/4104, -11.0/40}
};

static double rkf45_mu[RKF45_STAGE] = {
  16.0/135, 0.0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55
};

static double rkf45_mu2[RKF45_STAGE] = {
  25.0/216, 0.0, 1408.0/2565, 2197.0/4104, -1.0/5
};

double sum(int n, double *x)
{
  int i;
  double s = 0;
  for (i = 0; i < n; i++) s += x[i];
  return s;
}

/* $B4JC1$J%A%'%C%/(B */
void check_table()
{
  int i;
  for (i = 0; i < RKF45_STAGE; i++)
    printf("%g %g\n", rkf45_alpha[i], sum(i, rkf45_beta[i]));
  printf("rkf45_mu$B$NOB(B=%g\n", sum(RKF45_STAGE, rkf45_mu));
  printf("rkf45_mu2$B$NOB(B=%g\n", sum(RKF45_STAGE, rkf45_mu2));
}

double dotprod(int n, double *x, double *y)
{
  int i;
  double s = 0;
  for (i = 0; i < n; i++) s += x[i] * y[i];
  return s;
}

/*
 * RKF45 $B$G;~9o(B *t $B$+$i!"(B1$B%9%F%C%W?J$`(B
 * $BFC$KLdBj$J$1$l$P(B *hh $B$@$1?J$`$,!"C10LD9$5$"$?$j$N5vMF8m:98B3&(B eps_tol
 * ($B2r@b$N&E(B/H$B$KAjEv(B) $B$rC#@.$9$k$?$a$K!"I,MW$J$i$P9o$_I}$N<+F0D4@a$r$9$k!#(B
 */
int rkf45(double *newx,
	  double *t, double x, double *hh, double hmin, double eps_tol,
	  double f(double, double))
{
  int i, s = RKF45_STAGE;
  double k[RKF45_STAGE], newx1, newx2, h = *hh, error;
  do {
    /* k1,k2,k3,k4,k5,k6 $B$r7W;;$9$k(B */
    for (i = 0; i < s; i++)
      k[i] = h * f(*t + rkf45_alpha[i] * h, x + dotprod(i, rkf45_beta[i], k));
    /* 5$B.$5$/$J$i$J$$(B\n");
      exit(1);
    }
    printf("t=%25.16f, x=%25.16g\n", t, x);
  }
  printf("%25.16g\n", x);
  return 0;
}

$B$3$N%W%m%0%i%`$G$OGzH/$NNc$K=P$7$?=i4|CMLdBj(B

$$\frac{\Dx}{\Dt}=x^2,\quad x(0)=1$$

$B$r2r$$$F$$$k$,!"(B ($B¹²) $BGzH/;~9o(B $t=1$ $B$N$H$3$m$G!"9o$_I}$,$I$s$I$s>.$5$/$J$C$F$$$/$N$,J,$+$k!#(B

$B$J$*!"6qBNE*$J;~4V9o$_I}$N@)8f$K$D$$$F6qBNE*$K=q$$$F$"$kK\$OB?$/$J$$(B ($B$=$N$?$a$b$"$C$F!"(B $B;d<+?H$h$/J,$+$C$F$$$J$$ItJ,$b$"$j!">e5-$N%W%m%0%i%`[email protected]$G$O$J$$(B)$B!#(B $B$=$NE@$G?9(B [4] $B$O5.=E$G$"$k!#(B $B$=$3$K7G:\$5$l$F$$$k%W%m%0%i%`$N%"%k%4%j%:%`$O(B (FORTRAN $B%W%m%0%i%`$r(B C $BIw$K=q$-D>$7$F<($;$P(B) $B0J2<$N$h$&$J$b$N$G$"$k(B ($B$H;W$&!#FI$_0c$$$,L5$1$l$P!#(B)$B!#(B $B$J$k$[$I$H;W$&$H$3$m$b$"$k$,!":#0l$DG

  maxiter = 100;
  epsmac = 1e-15;
  // t0 $B$+$i(B tn $B$^$G2r$/(B
  // eps $B$O5vMF@dBP8m:9!"(Bepsmac $B$O7W;;5!%$%W%7%m%sDxEY$N?t(B
  epsv = max(eps, epsmac); // $B$O$F!"AjBP8m:9$J$i$H$b$+$/!"@dBP8m:9$G6Z$,DL$k(B?
  // eb $B$OC10L;~4VEv$j5v$5$l$k8m:9(B
  eb = fabs(epsv / (tn - t0));
  //
  t = t0;
  h = tn - t0;
  //
  for (iter = 1; iter <= itermax; iter++) {
    // $B$H$K$+$/9o$_I}(B h $B$G?J$s$G$_$F!"6I=jBG$A@Z$j8m:9$N?dDj$b$9$k(B
    rkf(t, h, x0, &xn, &errmax);
    // $B?dDj8m:9$,;~4V9o$_I}(B h $B$"$?$j$N5vMF8m:9$h$j$b>.$5$$$+$I$&$+%A%'%C%/(B
    if (errmax <= fabs(eb*h)) { // $B>.$5$$>l9g(B
      // $B;~9o$r?J$a$k(B
      t += h;
      // $B?J$s$@J,$@$1;D$j;~4V$O8:$i$9(B ($B:G=*;~9o(B - $B8=:_;~9o(B)
      h = tn - t;
      // $B;D$j;~4V$,$J$1$l$P(B ($B7W;;$7$-$C$?$i(B) $BLa$k(B
      if (fabs(h) <= epsmac)
        return;
      // $B$^$@;D$j;~4V$,$"$l$PB39T(B
      x0 = xn;
    }
    else {
      // $B?dDj8m:9$,;~4V9o$_I}(B h $B$"$?$j$N5vMF8m:9$h$j$bBg$-$$>l9g(B
      // $B?dDj8m:9$H5vMF@dBP8m:9$NBg>.$rHf3S(B
      if (errmax > epsv)
        // $B?dDj8m:9$,5vMF@dBP8m:9$h$jBg$-$1$l$P9o$_I}$r0l5$$K(B 10% $B$K=L$a$k(B
        h *= 0.1;
      else {
        h *= 0.9 * sqrt(sqrt(eb * fabs(h) / errmax));
        if (fabs(h) < epsmax)
          break;
        }
      }
    }
  }
  fprintf(stderr,
         "rkf45(): t=%g $B$G:$Fq$K=P2q$$$^$7$?!#7W;;$rCfCG$7$^$9!#(B", t);
  return 1;

($B$J$*!"$3$N5] $B$K$b2?$+=q$$$F$"$C$?$h$&$J5$$,$9$k!#(B $BDDi>9d$5$s$NGzH/$N%7%_%e%l!<%7%g%s$O$I$&$@$C$?$+$J!D(B $B$=$N$&$A$K0lEY$^$8$a$KD4$Y$F$_$h$&!#$^$"!"$3$l$OFH$j8@$_$?$$$J$b$N$G$9!#(B)