New ODE methods: Euler, Improved Euler, Adams-Bashford-Moulton

New ODE methods: Euler, Improved Euler, Adams-Bashford-Moulton

[alfonso_acosta_mail]

[-- Attachment #1: Type: text/plain, Size: 614 bytes --]

Hi:

We needed those methods for a college project and we implemented them 
with GSL, I know Euler and Improved Euler are merely didactic methods 
but maybe some students (like us) could avoid reinventing the wheel.

On the other hand, 4 step Adams-Bashford-Moulton method could be  usable 
in a real application. We used Runge-Kutta to obtain the remainig three 
initial steps.

I attached the 3 files just in case you want to include them in the GSL 
distribution

Cheers

Alfonso Acosta

PS: The three methods were tested without step control, so we are not 
sure wether the error we provide is correct or not.

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: abm.c --]
[-- Type: text/x-csrc; name="abm.c", Size: 8498 bytes --]

/* abm.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * Copyright (C) 2004 Daniel Rodríguez
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Adams-Bashford-Moulton */

/* Authors: Daniel Rodríguez
 */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

#define FALSE 0
#define TRUE !FALSE

typedef struct
{
  double valueT;
  double * valueY;
  double * valueF;
} node_t;

typedef struct
{
  node_t * TY;
  int inicializado;
} abm_state_t;

static void *
abm_alloc (size_t dim)
{
  abm_state_t *state = (abm_state_t *) malloc (sizeof (abm_state_t));

  if (state == NULL) {
    GSL_ERROR_NULL ("failed to allocate space for abm_state", GSL_ENOMEM);
  }
  state->inicializado = 3;

  state->TY = (node_t*) malloc (4 * sizeof (node_t));
  if (state->TY == NULL) {
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[0].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[0].valueY == NULL) {
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[1].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[1].valueY == NULL) {
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[2].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[2].valueY == NULL) {
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[3].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[3].valueY == NULL) {
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }

  state->TY[0].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[0].valueF == NULL) {
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[1].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[1].valueF == NULL) {
    free(state->TY[0].valueF);
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[2].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[2].valueF == NULL) {
    free(state->TY[1].valueF);
    free(state->TY[0].valueF);
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[3].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[3].valueF == NULL) {
    free(state->TY[2].valueF);
    free(state->TY[1].valueF);
    free(state->TY[0].valueF);
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }

  return state;
}


static int
abm_apply (void *vstate,
           size_t dim,
           double t,
           double h,
           double y[],
           double yerr[],
           const double dydt_in[],
           double dydt_out[], 
           const gsl_odeiv_system * sys)
{
  abm_state_t *state = (abm_state_t *) vstate;

  int * inicializado = &(state->inicializado);

  double temp, proxt;
  size_t i;
  int status = 0;
  int s;

  double * p = (double*) malloc(dim * sizeof(double));
  double * Ftmp = (double*) malloc(dim * sizeof(double));
  node_t *TY = state->TY;

  if (*inicializado) { /* Estamos en el primer paso, hacemos RK4 manualmente */
    TY[3-*inicializado].valueT = t;
    DBL_MEMCPY (TY[3-*inicializado].valueY, y, dim);
    s = GSL_ODEIV_FN_EVAL
      (sys, t, TY[3-*inicializado].valueY, TY[3-*inicializado].valueF);
    GSL_STATUS_UPDATE (&status, s);
    const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4;
    gsl_odeiv_step * step = gsl_odeiv_step_alloc(T, dim);
    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc(dim);
    int s2 = gsl_odeiv_evolve_apply(e, NULL, step, sys, &t, t+4*h, &h, y);
    GSL_STATUS_UPDATE (&status, s2);
    DBL_MEMCPY(TY[3-*inicializado+1].valueY, y, dim);
    TY[3-*inicializado+1].valueT = t;
    s = GSL_ODEIV_FN_EVAL
      (sys, t, TY[3-*inicializado+1].valueY, TY[3-*inicializado+1].valueF);
    GSL_STATUS_UPDATE (&status, s);
    for (i=0; i<dim; i++) {
      yerr[i] = h * TY[3-*inicializado+1].valueF[i];
      if (dydt_out != NULL)
	dydt_out[i] = TY[3-*inicializado+1].valueF[i];
    }
    (*inicializado)--;
  } else {
    for (i=0; i<dim; i++) {
      p[i] =
	TY[3].valueY[i] + (h/24)*(-9*TY[0].valueF[i]+37*TY[1].valueF[i]-
				  59*TY[2].valueF[i]+55*TY[3].valueF[i]);
    }
    proxt = TY[3].valueT + h;
    GSL_ODEIV_FN_EVAL (sys, proxt, p, Ftmp);
    GSL_STATUS_UPDATE (&status, s);
    for (i=0; i<dim; i++) {
      temp = (h/24)*(TY[1].valueF[i]-5*TY[2].valueF[i]+
		     19*TY[3].valueF[i]+9*Ftmp[i]);
      p[i] = TY[3].valueY[i] + temp;
      yerr[i] = h * temp;
      if (dydt_out != NULL)
	dydt_out[i] = temp;
    }
    DBL_MEMCPY(TY[0].valueY, TY[1].valueY, dim);
    DBL_MEMCPY(TY[1].valueY, TY[2].valueY, dim);
    DBL_MEMCPY(TY[2].valueY, TY[3].valueY, dim);
    DBL_MEMCPY(TY[3].valueY, p, dim);
    TY[0].valueT = TY[1].valueT;
    TY[1].valueT = TY[2].valueT;
    TY[2].valueT = TY[3].valueT;
    TY[3].valueT = proxt;
    DBL_MEMCPY(TY[0].valueF, TY[1].valueF, dim);
    DBL_MEMCPY(TY[1].valueF, TY[2].valueF, dim);
    DBL_MEMCPY(TY[2].valueF, TY[3].valueF, dim);
    GSL_ODEIV_FN_EVAL (sys, proxt, TY[3].valueY, Ftmp);
    GSL_STATUS_UPDATE (&status, s);
    DBL_MEMCPY(TY[3].valueF, Ftmp, dim);
    
    DBL_MEMCPY(y, TY[3].valueY, dim);
  }
  return status;
}

static int
abm_reset (void *vstate, size_t dim)
{
  abm_state_t *state = (abm_state_t *) vstate;

  DBL_ZERO_MEMSET (state->TY[0].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[1].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[2].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[3].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[0].valueF, dim);
  DBL_ZERO_MEMSET (state->TY[1].valueF, dim);
  DBL_ZERO_MEMSET (state->TY[2].valueF, dim);
  DBL_ZERO_MEMSET (state->TY[3].valueF, dim);
  state->TY[0].valueT = 0;
  state->TY[1].valueT = 0;
  state->TY[2].valueT = 0;
  state->TY[3].valueT = 0;
  state->inicializado = 3;

  return GSL_SUCCESS;
}

static unsigned int
abm_order (void *vstate)
{
  abm_state_t *state = (abm_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 4;
}

static void
abm_free (void *vstate)
{
  abm_state_t *state = (abm_state_t *) vstate;
  free (state->TY[0].valueY);
  free (state->TY[1].valueY);
  free (state->TY[2].valueY);
  free (state->TY[3].valueY);
  free (state->TY[0].valueF);
  free (state->TY[1].valueF);
  free (state->TY[2].valueF);
  free (state->TY[3].valueF);
  free (state->TY);
  free (state);
}

static const gsl_odeiv_step_type abm_type = { "abm",    /* name */
  0,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &abm_alloc,
  &abm_apply,
  &abm_reset,
  &abm_order,
  &abm_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_abm = &abm_type;

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: euler.c --]
[-- Type: text/x-csrc; name="euler.c", Size: 3204 bytes --]

/* euler.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * Copyright (C) 2004 Alfonso Acosta & Daniel Rodríguez
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Euler Method*/

/* Authors:  Alfonso Acosta & Daniel Rodríguez
 */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k;
}
euler_state_t;

static void *
euler_alloc (size_t dim)
{
  euler_state_t *state = (euler_state_t *) malloc (sizeof (euler_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for euler_state", GSL_ENOMEM);
    }

  state->k = (double *) malloc (dim * sizeof (double));

  if (state->k == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
    }

  return state;
}


static int
euler_apply (void *vstate,
           size_t dim,
           double t,
           double h,
           double y[],
           double yerr[],
           const double dydt_in[],
           double dydt_out[], 
           const gsl_odeiv_system * sys)
{
  euler_state_t *state = (euler_state_t *) vstate;
  
  double temp;
  size_t i;
  int status = 0;

  double *const k = state->k;


  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k);
      GSL_STATUS_UPDATE (&status, s);
    }
  

  for (i = 0; i < dim; i++)
    {
      temp = y[i]; /* save y[i] */ 
      y[i] = h * k[i];   
      yerr[i] = h * y[i];
      y[i] += temp;
      if (dydt_out != NULL)
          dydt_out[i] = k[i];
    }

  return status;
}

static int
euler_reset (void *vstate, size_t dim)
{
  euler_state_t *state = (euler_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k, dim);

  return GSL_SUCCESS;
}

static unsigned int
euler_order (void *vstate)
{
  euler_state_t *state = (euler_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 1;
}

static void
euler_free (void *vstate)
{
  euler_state_t *state = (euler_state_t *) vstate;
  free (state->k);
  free (state);
}

static const gsl_odeiv_step_type euler_type = { "euler",    /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &euler_alloc,
  &euler_apply,
  &euler_reset,
  &euler_order,
  &euler_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_euler = &euler_type;

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #4: eulerplus.c --]
[-- Type: text/x-csrc; name="eulerplus.c", Size: 4260 bytes --]

/* eulerplus.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * Copyright (C) 2004 Alfonso Acosta & Daniel Rodríguez
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Improved Euler Method */

/* Authors:  Alfonso Acosta & Daniel Rodríguez
 */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k1;
  double *k2;
  double *ytmp;
}
eulerplus_state_t;

static void *
eulerplus_alloc (size_t dim)
{
  eulerplus_state_t *state =
    (eulerplus_state_t *) malloc (sizeof (eulerplus_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for eulerplus_state",
		      GSL_ENOMEM);
    }

  state->k1 = (double *) malloc (dim * sizeof (double));

  if (state->k1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
    }

  state->k2 = (double *) malloc (dim * sizeof (double));

  if (state->k2 == 0)
    {
      free (state->k1);    
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
    }

  
  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  return state;
}


static int
eulerplus_apply (void *vstate,
           size_t dim,
           double t,
           double h,
           double y[],
           double yerr[],
           const double dydt_in[],
           double dydt_out[], 
           const gsl_odeiv_system * sys)
{
  eulerplus_state_t *state = (eulerplus_state_t *) vstate;

  double temp;
  size_t i;
  int status = 0;

  double *const k1 = state->k1;
  double *const k2 = state->k2;
  double *const ytmp = state->ytmp;


  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k1, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
      GSL_STATUS_UPDATE (&status, s);
    }
  /* aplicamos euler */
  for (i = 0; i < dim; i++)
    {
      ytmp[i] = y[i] + h * k1[i] ;
    }
  
  {
      int s = GSL_ODEIV_FN_EVAL (sys, t+h, ytmp, k2);
      GSL_STATUS_UPDATE (&status, s);	       
  }

  for (i = 0; i < dim; i++)
   {  
      temp = (k1[i]+k2[i]) / 2; /* la pendiente en x + h */ 
      if (dydt_out != NULL)
          dydt_out[i] = temp;
      yerr[i] = h * temp;
      y[i] += h*temp;
    }


  return status;
}

static int
eulerplus_reset (void *vstate, size_t dim)
{
  eulerplus_state_t *state = (eulerplus_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k1, dim);
  DBL_ZERO_MEMSET (state->k2, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);

  return GSL_SUCCESS;
}

static unsigned int
eulerplus_order (void *vstate)
{
  eulerplus_state_t *state = (eulerplus_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 2;
}

static void
eulerplus_free (void *vstate)
{
  eulerplus_state_t *state = (eulerplus_state_t *) vstate;
  free (state->k1);
  free (state->k2);
  free (state->ytmp);
  free (state);
}

static const gsl_odeiv_step_type eulerplus_type = { "eulerplus",    /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &eulerplus_alloc,
  &eulerplus_apply,
  &eulerplus_reset,
  &eulerplus_order,
  &eulerplus_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_eulerplus = &eulerplus_type;

Re: New ODE methods: Euler, Improved Euler, Adams-Bashford-Moulton

[alfonso_acosta_mail]

[-- Attachment #1: Type: text/plain, Size: 66 bytes --]

Sorry, I forgot to correct some comments in spanish from abm.c




[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: abm.c --]
[-- Type: text/x-csrc; name="abm.c", Size: 8473 bytes --]

/* abm.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * Copyright (C) 2004 Daniel Rodríguez
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* 4 step Adams-Bashford-Moulton */

/* Authors: Daniel Rodríguez
 */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

#define FALSE 0
#define TRUE !FALSE

typedef struct
{
  double valueT;
  double * valueY;
  double * valueF;
} node_t;

typedef struct
{
  node_t * TY;
  int initialized;
} abm_state_t;

static void *
abm_alloc (size_t dim)
{
  abm_state_t *state = (abm_state_t *) malloc (sizeof (abm_state_t));

  if (state == NULL) {
    GSL_ERROR_NULL ("failed to allocate space for abm_state", GSL_ENOMEM);
  }
  state->initialized = 3;

  state->TY = (node_t*) malloc (4 * sizeof (node_t));
  if (state->TY == NULL) {
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[0].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[0].valueY == NULL) {
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[1].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[1].valueY == NULL) {
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[2].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[2].valueY == NULL) {
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[3].valueY = (double*) malloc (dim * sizeof (double));
  if (state->TY[3].valueY == NULL) {
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }

  state->TY[0].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[0].valueF == NULL) {
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[1].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[1].valueF == NULL) {
    free(state->TY[0].valueF);
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[2].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[2].valueF == NULL) {
    free(state->TY[1].valueF);
    free(state->TY[0].valueF);
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }
  state->TY[3].valueF = (double*) malloc (dim * sizeof (double));
  if (state->TY[3].valueF == NULL) {
    free(state->TY[2].valueF);
    free(state->TY[1].valueF);
    free(state->TY[0].valueF);
    free(state->TY[3].valueY);
    free(state->TY[2].valueY);
    free(state->TY[1].valueY);
    free(state->TY[0].valueY);
    free(state);
    GSL_ERROR_NULL("failed to allocate space for valueY", GSL_ENOMEM);
  }

  return state;
}


static int
abm_apply (void *vstate,
           size_t dim,
           double t,
           double h,
           double y[],
           double yerr[],
           const double dydt_in[],
           double dydt_out[], 
           const gsl_odeiv_system * sys)
{
  abm_state_t *state = (abm_state_t *) vstate;

  int * initialized = &(state->initialized);

  double temp, proxt;
  size_t i;
  int status = 0;
  int s;

  double * p = (double*) malloc(dim * sizeof(double));
  double * Ftmp = (double*) malloc(dim * sizeof(double));
  node_t *TY = state->TY;

  if (*initialized) { /* obtain the first 3 steps with RK4 */
    TY[3-*initialized].valueT = t;
    DBL_MEMCPY (TY[3-*initialized].valueY, y, dim);
    s = GSL_ODEIV_FN_EVAL
      (sys, t, TY[3-*initialized].valueY, TY[3-*initialized].valueF);
    GSL_STATUS_UPDATE (&status, s);
    const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4;
    gsl_odeiv_step * step = gsl_odeiv_step_alloc(T, dim);
    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc(dim);
    int s2 = gsl_odeiv_evolve_apply(e, NULL, step, sys, &t, t+4*h, &h, y);
    GSL_STATUS_UPDATE (&status, s2);
    DBL_MEMCPY(TY[3-*initialized+1].valueY, y, dim);
    TY[3-*initialized+1].valueT = t;
    s = GSL_ODEIV_FN_EVAL
      (sys, t, TY[3-*initialized+1].valueY, TY[3-*initialized+1].valueF);
    GSL_STATUS_UPDATE (&status, s);
    for (i=0; i<dim; i++) {
      yerr[i] = h * TY[3-*initialized+1].valueF[i];
      if (dydt_out != NULL)
	dydt_out[i] = TY[3-*initialized+1].valueF[i];
    }
    (*initialized)--;
  } else {
    for (i=0; i<dim; i++) {
      p[i] =
	TY[3].valueY[i] + (h/24)*(-9*TY[0].valueF[i]+37*TY[1].valueF[i]-
				  59*TY[2].valueF[i]+55*TY[3].valueF[i]);
    }
    proxt = TY[3].valueT + h;
    GSL_ODEIV_FN_EVAL (sys, proxt, p, Ftmp);
    GSL_STATUS_UPDATE (&status, s);
    for (i=0; i<dim; i++) {
      temp = (h/24)*(TY[1].valueF[i]-5*TY[2].valueF[i]+
		     19*TY[3].valueF[i]+9*Ftmp[i]);
      p[i] = TY[3].valueY[i] + temp;
      yerr[i] = h * temp;
      if (dydt_out != NULL)
	dydt_out[i] = temp;
    }
    DBL_MEMCPY(TY[0].valueY, TY[1].valueY, dim);
    DBL_MEMCPY(TY[1].valueY, TY[2].valueY, dim);
    DBL_MEMCPY(TY[2].valueY, TY[3].valueY, dim);
    DBL_MEMCPY(TY[3].valueY, p, dim);
    TY[0].valueT = TY[1].valueT;
    TY[1].valueT = TY[2].valueT;
    TY[2].valueT = TY[3].valueT;
    TY[3].valueT = proxt;
    DBL_MEMCPY(TY[0].valueF, TY[1].valueF, dim);
    DBL_MEMCPY(TY[1].valueF, TY[2].valueF, dim);
    DBL_MEMCPY(TY[2].valueF, TY[3].valueF, dim);
    GSL_ODEIV_FN_EVAL (sys, proxt, TY[3].valueY, Ftmp);
    GSL_STATUS_UPDATE (&status, s);
    DBL_MEMCPY(TY[3].valueF, Ftmp, dim);
    
    DBL_MEMCPY(y, TY[3].valueY, dim);
  }
  return status;
}

static int
abm_reset (void *vstate, size_t dim)
{
  abm_state_t *state = (abm_state_t *) vstate;

  DBL_ZERO_MEMSET (state->TY[0].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[1].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[2].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[3].valueY, dim);
  DBL_ZERO_MEMSET (state->TY[0].valueF, dim);
  DBL_ZERO_MEMSET (state->TY[1].valueF, dim);
  DBL_ZERO_MEMSET (state->TY[2].valueF, dim);
  DBL_ZERO_MEMSET (state->TY[3].valueF, dim);
  state->TY[0].valueT = 0;
  state->TY[1].valueT = 0;
  state->TY[2].valueT = 0;
  state->TY[3].valueT = 0;
  state->initialized = 3;

  return GSL_SUCCESS;
}

static unsigned int
abm_order (void *vstate)
{
  abm_state_t *state = (abm_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 4;
}

static void
abm_free (void *vstate)
{
  abm_state_t *state = (abm_state_t *) vstate;
  free (state->TY[0].valueY);
  free (state->TY[1].valueY);
  free (state->TY[2].valueY);
  free (state->TY[3].valueY);
  free (state->TY[0].valueF);
  free (state->TY[1].valueF);
  free (state->TY[2].valueF);
  free (state->TY[3].valueF);
  free (state->TY);
  free (state);
}

static const gsl_odeiv_step_type abm_type = { "abm",    /* name */
  0,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &abm_alloc,
  &abm_apply,
  &abm_reset,
  &abm_order,
  &abm_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_abm = &abm_type;

Re: New ODE methods: Euler, Improved Euler, Adams-Bashford-Moulton

[Brian Gough]

alfonso_acosta_mail@yahoo.es writes:
 > We needed those methods for a college project and we implemented them 
 > with GSL, I know Euler and Improved Euler are merely didactic methods 
 > but maybe some students (like us) could avoid reinventing the wheel.
 > 
 > On the other hand, 4 step Adams-Bashford-Moulton method could be  usable 
 > in a real application. We used Runge-Kutta to obtain the remainig three 
 > initial steps.
 > 
 > I attached the 3 files just in case you want to include them in the GSL 
 > distribution

Thanks, I will take a look at them.

-- 
Brian Gough

gsl_root_test_residual(double f, double epsabs)

[Karsten Howes]

Hi,

I am working with the 1D root finder framework and I am trying to figure 
out how to get the value of my function at the current best root estimate 
(without re-evaluating my function of course).  I need it because I want 
to stop when the function is sufficiently close to zero and I dont 
actually care about the exact location of the root.  The information does 
not appear to be in the solver state:

typedef struct
   {
     const gsl_root_fsolver_type * type;
     gsl_function * function ;
     double root ;
     double x_lower;
     double x_upper;
     void *state;
   }
gsl_root_fsolver;

How do I use the function gsl_root_test_residual?

If it doesn't work, one possibility is for me to maintain a map x->f for 
all previous valuations of f and then just look up f in the map using the 
current best root.  Sorry if I'm missing something obvious.

Thanks for a great library Brian et al.

Best regards,
Karsten

Re: gsl_root_test_residual(double f, double epsabs)

[Brian Gough]

Karsten Howes writes:
 > I am working with the 1D root finder framework and I am trying to figure 
 > out how to get the value of my function at the current best root estimate 
 > (without re-evaluating my function of course).  

You do need to evaluate the function f(root).  Depending on the
algorithm it may not be a reevaluation (see e.g. bisection).

-- 
Brian Gough

Network Theory Ltd,
Commercial support for GSL --- http://www.network-theory.co.uk/gsl/

Re: gsl_root_test_residual(double f, double epsabs)

[Karsten Howes]

OK, thanks Brian.  I am using Brent's method which appears always to have 
valued the function at the root.  I solved the problem by storing a map 
root->value inside my function.  I have included the code below.  It's a 
C++/stl/boost wrapper for the brent solver, which will be of no interest 
to those who are hard core C programmers. And, it has a few pointers of 
overhead, etc, which doesnt matter to me in my applications.  It does have 
the advantage of being able to find roots of member functions though 
(using e.g. boost::bind).

Best regards,
Karsten

   //params for gsl root finder.  Includes the function
   struct DoubleGslBoostFuncWrapper {
     boost::function<double (double)> f;
     std::map<double, double> f_map;
   };

   static double doubleGslFunc(double x, void *boostFunc){
     double f = ((DoubleGslBoostFuncWrapper *) boostFunc)->f(x);
     ((DoubleGslBoostFuncWrapper *) boostFunc)->f_map[x] = f;
     return f;
   }

   //Finds root of f with error of tolerance
   int brent(boost::function<double (double)> &f,
	    double x_lo, double x_hi, double tolerance, double &root){
     int status;
     int iter = 0, max_iter = 20;
     const gsl_root_fsolver_type *T;
     gsl_root_fsolver *s;
     gsl_function F;

     DoubleGslBoostFuncWrapper gslBoostFuncWrapper;
     gslBoostFuncWrapper.f = f;
     F.function = &doubleGslFunc;
     F.params = &gslBoostFuncWrapper;

     T = gsl_root_fsolver_brent;
     s = gsl_root_fsolver_alloc (T);
     gsl_root_fsolver_set (s, &F, x_lo, x_hi);

     printf ("using %s method\n",
	    gsl_root_fsolver_name (s));

     printf ("%5s [%9s, %9s] %9s %10s %9s\n",
	    "iter", "lower", "upper", "root",
	    "err", "err(est)");

     do
       {
	iter++;
	status = gsl_root_fsolver_iterate (s);
	root = gsl_root_fsolver_root (s);
	x_lo = gsl_root_fsolver_x_lower (s);
	x_hi = gsl_root_fsolver_x_upper (s);
	//status = gsl_root_test_interval (x_lo, x_hi, 0, 0.001);
	double residual = gslBoostFuncWrapper.f_map[root];
	status = gsl_root_test_residual (residual, tolerance);

	if (status == GSL_SUCCESS)
	  printf ("Converged:\n");

	printf ("%5d [%.7f, %.7f] %.7f %.7f, f=%.7f\n",
		iter, x_lo, x_hi, root, x_hi - x_lo, residual);
       }
     while (status == GSL_CONTINUE && iter < max_iter);
     return status;
   }



On Wed, 08 Sep 2004 18:26:59 +0100, Brian Gough <bjg@network-theory.co.uk> 
wrote:

> Karsten Howes writes:
>  > I am working with the 1D root finder framework and I am trying to 
> figure
>  > out how to get the value of my function at the current best root 
> estimate
>  > (without re-evaluating my function of course).
>
> You do need to evaluate the function f(root).  Depending on the
> algorithm it may not be a reevaluation (see e.g. bisection).
>