Annotation of OpenXM/src/hgm/gsl-t-1/src/t-cstd.c, Revision 1.1
1.1 ! takayama 1: /* ode-initval/cstd.c
! 2: *
! 3: * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
! 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: #include "config.h"
! 21: #include <stdlib.h>
! 22: #include <math.h>
! 23: #include "t-gsl_errno.h"
! 24: #include "t-gsl_math.h"
! 25: #include "t-gsl_odeiv.h"
! 26:
! 27: typedef struct
! 28: {
! 29: double eps_abs;
! 30: double eps_rel;
! 31: double a_y;
! 32: double a_dydt;
! 33: }
! 34: std_control_state_t;
! 35:
! 36: static void *
! 37: std_control_alloc (void)
! 38: {
! 39: std_control_state_t * s =
! 40: (std_control_state_t *) malloc (sizeof(std_control_state_t));
! 41:
! 42: if (s == 0)
! 43: {
! 44: GSL_ERROR_NULL ("failed to allocate space for std_control_state",
! 45: GSL_ENOMEM);
! 46: }
! 47:
! 48: return s;
! 49: }
! 50:
! 51: static int
! 52: std_control_init (void * vstate,
! 53: double eps_abs, double eps_rel, double a_y, double a_dydt)
! 54: {
! 55: std_control_state_t * s = (std_control_state_t *) vstate;
! 56:
! 57: if (eps_abs < 0)
! 58: {
! 59: GSL_ERROR ("eps_abs is negative", GSL_EINVAL);
! 60: }
! 61: else if (eps_rel < 0)
! 62: {
! 63: GSL_ERROR ("eps_rel is negative", GSL_EINVAL);
! 64: }
! 65: else if (a_y < 0)
! 66: {
! 67: GSL_ERROR ("a_y is negative", GSL_EINVAL);
! 68: }
! 69: else if (a_dydt < 0)
! 70: {
! 71: GSL_ERROR ("a_dydt is negative", GSL_EINVAL);
! 72: }
! 73:
! 74: s->eps_rel = eps_rel;
! 75: s->eps_abs = eps_abs;
! 76: s->a_y = a_y;
! 77: s->a_dydt = a_dydt;
! 78:
! 79: return GSL_SUCCESS;
! 80: }
! 81:
! 82: static int
! 83: std_control_hadjust(void * vstate, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h)
! 84: {
! 85: std_control_state_t *state = (std_control_state_t *) vstate;
! 86:
! 87: const double eps_abs = state->eps_abs;
! 88: const double eps_rel = state->eps_rel;
! 89: const double a_y = state->a_y;
! 90: const double a_dydt = state->a_dydt;
! 91:
! 92: const double S = 0.9;
! 93: const double h_old = *h;
! 94:
! 95: double rmax = DBL_MIN;
! 96: size_t i;
! 97:
! 98: for(i=0; i<dim; i++) {
! 99: const double D0 =
! 100: eps_rel * (a_y * fabs(y[i]) + a_dydt * fabs(h_old * yp[i])) + eps_abs;
! 101: const double r = fabs(yerr[i]) / fabs(D0);
! 102: // rmax = GSL_MAX_DBL(r, rmax);
! 103: rmax = (r > rmax? r: rmax);
! 104: }
! 105:
! 106: if(rmax > 1.1) {
! 107: /* decrease step, no more than factor of 5, but a fraction S more
! 108: than scaling suggests (for better accuracy) */
! 109: double r = S / pow(rmax, 1.0/ord);
! 110:
! 111: if (r < 0.2)
! 112: r = 0.2;
! 113:
! 114: *h = r * h_old;
! 115:
! 116: return GSL_ODEIV_HADJ_DEC;
! 117: }
! 118: else if(rmax < 0.5) {
! 119: /* increase step, no more than factor of 5 */
! 120: double r = S / pow(rmax, 1.0/(ord+1.0));
! 121:
! 122: if (r > 5.0)
! 123: r = 5.0;
! 124:
! 125: if (r < 1.0) /* don't allow any decrease caused by S<1 */
! 126: r = 1.0;
! 127:
! 128: *h = r * h_old;
! 129:
! 130: return GSL_ODEIV_HADJ_INC;
! 131: }
! 132: else {
! 133: /* no change */
! 134: return GSL_ODEIV_HADJ_NIL;
! 135: }
! 136: }
! 137:
! 138:
! 139: static void
! 140: std_control_free (void * vstate)
! 141: {
! 142: std_control_state_t *state = (std_control_state_t *) vstate;
! 143: free (state);
! 144: }
! 145:
! 146: static const gsl_odeiv_control_type std_control_type =
! 147: {"standard", /* name */
! 148: &std_control_alloc,
! 149: &std_control_init,
! 150: &std_control_hadjust,
! 151: &std_control_free};
! 152:
! 153: const gsl_odeiv_control_type *gsl_odeiv_control_standard = &std_control_type;
! 154:
! 155:
! 156: gsl_odeiv_control *
! 157: gsl_odeiv_control_standard_new(double eps_abs, double eps_rel,
! 158: double a_y, double a_dydt)
! 159: {
! 160: gsl_odeiv_control * c =
! 161: gsl_odeiv_control_alloc (gsl_odeiv_control_standard);
! 162:
! 163: int status = gsl_odeiv_control_init (c, eps_abs, eps_rel, a_y, a_dydt);
! 164:
! 165: if (status != GSL_SUCCESS)
! 166: {
! 167: gsl_odeiv_control_free (c);
! 168: GSL_ERROR_NULL ("error trying to initialize control", status);
! 169: }
! 170:
! 171: return c;
! 172: }
! 173:
! 174: gsl_odeiv_control *
! 175: gsl_odeiv_control_y_new(double eps_abs, double eps_rel)
! 176: {
! 177: return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 1.0, 0.0);
! 178: }
! 179:
! 180:
! 181: gsl_odeiv_control *
! 182: gsl_odeiv_control_yp_new(double eps_abs, double eps_rel)
! 183: {
! 184: return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 0.0, 1.0);
! 185: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to t-cstd.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