Annotation of OpenXM/src/hgm/gsl-t-1/src/t-rkf45.c, Revision 1.1
1.1 ! takayama 1: /* ode-initval/rkf45.c
! 2: *
! 3: * Copyright (C) 2001, 2004, 2007 Brian Gough
! 4: *
! 5: * This program is free software; you can redistribute it and/or modify
! 6: * it under the terms of the GNU General Public License as published by
! 7: * the Free Software Foundation; either version 3 of the License, or (at
! 8: * your option) any later version.
! 9: *
! 10: * This program is distributed in the hope that it will be useful, but
! 11: * WITHOUT ANY WARRANTY; without even the implied warranty of
! 12: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! 13: * General Public License for more details.
! 14: *
! 15: * You should have received a copy of the GNU General Public License
! 16: * along with this program; if not, write to the Free Software
! 17: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
! 18: */
! 19:
! 20: /* Runge-Kutta-Fehlberg 4(5)*/
! 21:
! 22: /* Reference eg. Hairer, E., Norsett S.P., Wanner, G. Solving ordinary
! 23: differential equations I, Nonstiff Problems, 2nd revised edition,
! 24: Springer, 2000.
! 25: */
! 26:
! 27: #include "config.h"
! 28: #include <stdlib.h>
! 29: #include <string.h>
! 30: #include "t-gsl_errno.h"
! 31: #include "t-gsl_odeiv.h"
! 32:
! 33: #include "odeiv_util.h"
! 34:
! 35: /* Runge-Kutta-Fehlberg coefficients. Zero elements left out */
! 36:
! 37: static const double ah[] = { 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 };
! 38: static const double b3[] = { 3.0/32.0, 9.0/32.0 };
! 39: static const double b4[] = { 1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0};
! 40: static const double b5[] = { 8341.0/4104.0, -32832.0/4104.0, 29440.0/4104.0, -845.0/4104.0};
! 41: static const double b6[] = { -6080.0/20520.0, 41040.0/20520.0, -28352.0/20520.0, 9295.0/20520.0, -5643.0/20520.0};
! 42:
! 43: static const double c1 = 902880.0/7618050.0;
! 44: static const double c3 = 3953664.0/7618050.0;
! 45: static const double c4 = 3855735.0/7618050.0;
! 46: static const double c5 = -1371249.0/7618050.0;
! 47: static const double c6 = 277020.0/7618050.0;
! 48:
! 49: /* These are the differences of fifth and fourth order coefficients
! 50: for error estimation */
! 51:
! 52: static const double ec[] = { 0.0,
! 53: 1.0 / 360.0,
! 54: 0.0,
! 55: -128.0 / 4275.0,
! 56: -2197.0 / 75240.0,
! 57: 1.0 / 50.0,
! 58: 2.0 / 55.0
! 59: };
! 60:
! 61: typedef struct
! 62: {
! 63: double *k1;
! 64: double *k2;
! 65: double *k3;
! 66: double *k4;
! 67: double *k5;
! 68: double *k6;
! 69: double *y0;
! 70: double *ytmp;
! 71: }
! 72: rkf45_state_t;
! 73:
! 74: static void *
! 75: rkf45_alloc (size_t dim)
! 76: {
! 77: rkf45_state_t *state = (rkf45_state_t *) malloc (sizeof (rkf45_state_t));
! 78:
! 79: if (state == 0)
! 80: {
! 81: GSL_ERROR_NULL ("failed to allocate space for rkf45_state", GSL_ENOMEM);
! 82: }
! 83:
! 84: state->k1 = (double *) malloc (dim * sizeof (double));
! 85:
! 86: if (state->k1 == 0)
! 87: {
! 88: free (state);
! 89: GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
! 90: }
! 91:
! 92: state->k2 = (double *) malloc (dim * sizeof (double));
! 93:
! 94: if (state->k2 == 0)
! 95: {
! 96: free (state->k1);
! 97: free (state);
! 98: GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
! 99: }
! 100:
! 101: state->k3 = (double *) malloc (dim * sizeof (double));
! 102:
! 103: if (state->k3 == 0)
! 104: {
! 105: free (state->k2);
! 106: free (state->k1);
! 107: free (state);
! 108: GSL_ERROR_NULL ("failed to allocate space for k3", GSL_ENOMEM);
! 109: }
! 110:
! 111: state->k4 = (double *) malloc (dim * sizeof (double));
! 112:
! 113: if (state->k4 == 0)
! 114: {
! 115: free (state->k3);
! 116: free (state->k2);
! 117: free (state->k1);
! 118: free (state);
! 119: GSL_ERROR_NULL ("failed to allocate space for k4", GSL_ENOMEM);
! 120: }
! 121:
! 122: state->k5 = (double *) malloc (dim * sizeof (double));
! 123:
! 124: if (state->k5 == 0)
! 125: {
! 126: free (state->k4);
! 127: free (state->k3);
! 128: free (state->k2);
! 129: free (state->k1);
! 130: free (state);
! 131: GSL_ERROR_NULL ("failed to allocate space for k5", GSL_ENOMEM);
! 132: }
! 133:
! 134: state->k6 = (double *) malloc (dim * sizeof (double));
! 135:
! 136: if (state->k6 == 0)
! 137: {
! 138: free (state->k5);
! 139: free (state->k4);
! 140: free (state->k3);
! 141: free (state->k2);
! 142: free (state->k1);
! 143: free (state);
! 144: GSL_ERROR_NULL ("failed to allocate space for k6", GSL_ENOMEM);
! 145: }
! 146:
! 147: state->y0 = (double *) malloc (dim * sizeof (double));
! 148:
! 149: if (state->y0 == 0)
! 150: {
! 151: free (state->k6);
! 152: free (state->k5);
! 153: free (state->k4);
! 154: free (state->k3);
! 155: free (state->k2);
! 156: free (state->k1);
! 157: free (state);
! 158: GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
! 159: }
! 160:
! 161: state->ytmp = (double *) malloc (dim * sizeof (double));
! 162:
! 163: if (state->ytmp == 0)
! 164: {
! 165: free (state->y0);
! 166: free (state->k6);
! 167: free (state->k5);
! 168: free (state->k4);
! 169: free (state->k3);
! 170: free (state->k2);
! 171: free (state->k1);
! 172: free (state);
! 173: GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
! 174: }
! 175:
! 176: return state;
! 177: }
! 178:
! 179:
! 180: static int
! 181: rkf45_apply (void *vstate,
! 182: size_t dim,
! 183: double t,
! 184: double h,
! 185: double y[],
! 186: double yerr[],
! 187: const double dydt_in[],
! 188: double dydt_out[], const gsl_odeiv_system * sys)
! 189: {
! 190: rkf45_state_t *state = (rkf45_state_t *) vstate;
! 191:
! 192: size_t i;
! 193:
! 194: double *const k1 = state->k1;
! 195: double *const k2 = state->k2;
! 196: double *const k3 = state->k3;
! 197: double *const k4 = state->k4;
! 198: double *const k5 = state->k5;
! 199: double *const k6 = state->k6;
! 200: double *const ytmp = state->ytmp;
! 201: double *const y0 = state->y0;
! 202:
! 203: DBL_MEMCPY (y0, y, dim);
! 204:
! 205: /* k1 step */
! 206: if (dydt_in != NULL)
! 207: {
! 208: DBL_MEMCPY (k1, dydt_in, dim);
! 209: }
! 210: else
! 211: {
! 212: int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
! 213:
! 214: if (s != GSL_SUCCESS)
! 215: {
! 216: return s;
! 217: }
! 218: }
! 219:
! 220: for (i = 0; i < dim; i++)
! 221: ytmp[i] = y[i] + ah[0] * h * k1[i];
! 222:
! 223: /* k2 step */
! 224: {
! 225: int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);
! 226:
! 227: if (s != GSL_SUCCESS)
! 228: {
! 229: return s;
! 230: }
! 231: }
! 232:
! 233: for (i = 0; i < dim; i++)
! 234: ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);
! 235:
! 236: /* k3 step */
! 237: {
! 238: int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);
! 239:
! 240: if (s != GSL_SUCCESS)
! 241: {
! 242: return s;
! 243: }
! 244: }
! 245:
! 246: for (i = 0; i < dim; i++)
! 247: ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[1] * k2[i] + b4[2] * k3[i]);
! 248:
! 249: /* k4 step */
! 250: {
! 251: int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);
! 252:
! 253: if (s != GSL_SUCCESS)
! 254: {
! 255: return s;
! 256: }
! 257: }
! 258:
! 259: for (i = 0; i < dim; i++)
! 260: ytmp[i] =
! 261: y[i] + h * (b5[0] * k1[i] + b5[1] * k2[i] + b5[2] * k3[i] +
! 262: b5[3] * k4[i]);
! 263:
! 264: /* k5 step */
! 265: {
! 266: int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);
! 267:
! 268: if (s != GSL_SUCCESS)
! 269: {
! 270: return s;
! 271: }
! 272: }
! 273:
! 274: for (i = 0; i < dim; i++)
! 275: ytmp[i] =
! 276: y[i] + h * (b6[0] * k1[i] + b6[1] * k2[i] + b6[2] * k3[i] +
! 277: b6[3] * k4[i] + b6[4] * k5[i]);
! 278:
! 279: /* k6 step and final sum */
! 280: {
! 281: int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);
! 282:
! 283: if (s != GSL_SUCCESS)
! 284: {
! 285: return s;
! 286: }
! 287: }
! 288:
! 289: for (i = 0; i < dim; i++)
! 290: {
! 291: const double d_i = c1 * k1[i] + c3 * k3[i] + c4 * k4[i] + c5 * k5[i] + c6 * k6[i];
! 292: y[i] += h * d_i;
! 293: }
! 294:
! 295: /* Derivatives at output */
! 296:
! 297: if (dydt_out != NULL)
! 298: {
! 299: int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);
! 300:
! 301: if (s != GSL_SUCCESS)
! 302: {
! 303: /* Restore initial values */
! 304: DBL_MEMCPY (y, y0, dim);
! 305:
! 306: return s;
! 307: }
! 308: }
! 309:
! 310: /* difference between 4th and 5th order */
! 311: for (i = 0; i < dim; i++)
! 312: {
! 313: yerr[i] = h * (ec[1] * k1[i] + ec[3] * k3[i] + ec[4] * k4[i]
! 314: + ec[5] * k5[i] + ec[6] * k6[i]);
! 315: }
! 316:
! 317: return GSL_SUCCESS;
! 318: }
! 319:
! 320:
! 321: static int
! 322: rkf45_reset (void *vstate, size_t dim)
! 323: {
! 324: rkf45_state_t *state = (rkf45_state_t *) vstate;
! 325:
! 326: DBL_ZERO_MEMSET (state->k1, dim);
! 327: DBL_ZERO_MEMSET (state->k2, dim);
! 328: DBL_ZERO_MEMSET (state->k3, dim);
! 329: DBL_ZERO_MEMSET (state->k4, dim);
! 330: DBL_ZERO_MEMSET (state->k5, dim);
! 331: DBL_ZERO_MEMSET (state->k6, dim);
! 332: DBL_ZERO_MEMSET (state->ytmp, dim);
! 333: DBL_ZERO_MEMSET (state->y0, dim);
! 334:
! 335: return GSL_SUCCESS;
! 336: }
! 337:
! 338: static unsigned int
! 339: rkf45_order (void *vstate)
! 340: {
! 341: rkf45_state_t *state = (rkf45_state_t *) vstate;
! 342: state = 0; /* prevent warnings about unused parameters */
! 343: return 5;
! 344: }
! 345:
! 346: static void
! 347: rkf45_free (void *vstate)
! 348: {
! 349: rkf45_state_t *state = (rkf45_state_t *) vstate;
! 350:
! 351: free (state->ytmp);
! 352: free (state->y0);
! 353: free (state->k6);
! 354: free (state->k5);
! 355: free (state->k4);
! 356: free (state->k3);
! 357: free (state->k2);
! 358: free (state->k1);
! 359: free (state);
! 360: }
! 361:
! 362: static const gsl_odeiv_step_type rkf45_type = { "rkf45", /* name */
! 363: 1, /* can use dydt_in */
! 364: 0, /* gives exact dydt_out */
! 365: &rkf45_alloc,
! 366: &rkf45_apply,
! 367: &rkf45_reset,
! 368: &rkf45_order,
! 369: &rkf45_free
! 370: };
! 371:
! 372: const gsl_odeiv_step_type *gsl_odeiv_step_rkf45 = &rkf45_type;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to t-rkf45.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / hgm / gsl-t-1 / src