Annotated Source Code

1/* ode-initval2/driver.c
2 * 
3 * Copyright (C) 2009, 2010 Tuomo Keskitalo
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/* Driver routine for odeiv2. This is a wrapper for low level GSL
21   functions that allows a simple interface to step, control and
22   evolve layers.
23 */
24
25#include <config.h>
26#include <math.h>
27#include <gsl/gsl_errno.h>
28#include <gsl/gsl_odeiv2.h>
29#include <gsl/gsl_machine.h>
30
31static gsl_odeiv2_driver *
32driver_alloc (const gsl_odeiv2_system * sys, const double hstart,
1Start Analysis.
33              const gsl_odeiv2_step_type * T)
34{
35  /* Allocates and initializes an ODE driver system. Step and evolve
36     objects are allocated here, but control object is allocated in
37     another function.
38   */
39
40  gsl_odeiv2_driver *state =
41    (gsl_odeiv2_driver *) malloc (sizeof (gsl_odeiv2_driver));
2Call a function.    malloc(sizeof(gsl_odeiv2_driver))
42
43  if (state == NULL)
3Take the false branch.    state == ((void *)0)
44    {
45      GSL_ERROR_NULL ("failed to allocate space for driver state",
46                      GSL_ENOMEM);
47    }
48
49  if (sys == NULL)
4Take the false branch.    sys == ((void *)0)
50    {
51      GSL_ERROR_NULL ("gsl_odeiv2_system must be defined", GSL_EINVAL);
52    }
53
54  {
55    const size_t dim = sys->dimension;
56
57    if (dim == 0)
5Take the false branch.    dim == 0
58      {
59        GSL_ERROR_NULL
60          ("gsl_odeiv2_system dimension must be a positive integer",
61           GSL_EINVAL);
62      }
63
64    state->sys = sys;
65
66    state->s = gsl_odeiv2_step_alloc (T, dim);
6Call a function.    gsl_odeiv2_step_alloc(T, dim)
67
68    if (state->s == NULL)
12Take the false branch.    state->s == ((void *)0)
69      {
70        free (state);
71        GSL_ERROR_NULL ("failed to allocate step object", GSL_ENOMEM);
72      }
73
74    state->e = gsl_odeiv2_evolve_alloc (dim);
13Call a function.    gsl_odeiv2_evolve_alloc(dim)
75  }
76
77  if (state->e == NULL)
26Take the false branch.    state->e == ((void *)0)
78    {
79      gsl_odeiv2_step_free (state->s);
80      free (state);
81      GSL_ERROR_NULL ("failed to allocate evolve object", GSL_ENOMEM);
82    }
83
84  if (hstart > 0.0 || hstart < 0.0)
27Take the false branch.    hstart > 0.
28Take the false branch.    hstart > 0. || hstart < 0.
85    {
86      state->h = hstart;
87    }
88  else
89    {
90      GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL);
29Call a function.    gsl_error("invalid hstart", "/home/xingming/workplace/experiment/fp/gsl-2.1/ode-initval2/driver.c", 90, GSL_EINVAL)
91    }
92
93  state->h = hstart;
94  state->hmin = 0.0;
95  state->hmax = GSL_DBL_MAX;
96  state->nmax = 0;
97  state->n = 0;
98  state->c = NULL;
99
100  return state;
101}
102
103int
104gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin)
105{
106  /* Sets minimum allowed step size fabs(hmin) for driver. It is
107     required that hmin <= fabs(h) <= hmax. */
108
109  if ((fabs (hmin) > fabs (d->h)) || (fabs (hmin) > d->hmax))
110    {
111      GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
112    }
113
114  d->hmin = fabs (hmin);
115
116  return GSL_SUCCESS;
117}
118
119int
120gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax)
121{
122  /* Sets maximum allowed step size fabs(hmax) for driver. It is
123     required that hmin <= fabs(h) <= hmax. */
124
125  if ((fabs (hmax) < fabs (d->h)) || (fabs (hmax) < d->hmin))
126    {
127      GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
128    }
129
130  if (hmax > 0.0 || hmax < 0.0)
131    {
132      d->hmax = fabs (hmax);
133    }
134  else
135    {
136      GSL_ERROR_NULL ("invalid hmax", GSL_EINVAL);
137    }
138
139  return GSL_SUCCESS;
140}
141
142int
143gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d,
144                            const unsigned long int nmax)
145{
146  /* Sets maximum number of allowed steps (nmax) for driver */
147
148  d->nmax = nmax;
149
150  return GSL_SUCCESS;
151}
152
153gsl_odeiv2_driver *
154gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * sys,
155                               const gsl_odeiv2_step_type * T,
156                               const double hstart,
157                               const double epsabs, const double epsrel)
158{
159  /* Initializes an ODE driver system with control object of type y_new. */
160
161  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
162
163  if (state == NULL)
164    {
165      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
166    }
167
168  if (epsabs >= 0.0 && epsrel >= 0.0)
169    {
170      state->c = gsl_odeiv2_control_y_new (epsabs, epsrel);
171
172      if (state->c == NULL)
173        {
174          gsl_odeiv2_driver_free (state);
175          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
176        }
177    }
178  else
179    {
180      gsl_odeiv2_driver_free (state);
181      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
182    }
183
184  /* Distribute pointer to driver object */
185
186  gsl_odeiv2_step_set_driver (state->s, state);
187  gsl_odeiv2_evolve_set_driver (state->e, state);
188  gsl_odeiv2_control_set_driver (state->c, state);
189
190  return state;
191}
192
193gsl_odeiv2_driver *
194gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * sys,
195                                const gsl_odeiv2_step_type * T,
196                                const double hstart,
197                                const double epsabs, const double epsrel)
198{
199  /* Initializes an ODE driver system with control object of type yp_new. */
200
201  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
202
203  if (state == NULL)
204    {
205      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
206    }
207
208  if (epsabs >= 0.0 && epsrel >= 0.0)
209    {
210      state->c = gsl_odeiv2_control_yp_new (epsabs, epsrel);
211
212      if (state->c == NULL)
213        {
214          gsl_odeiv2_driver_free (state);
215          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
216        }
217    }
218  else
219    {
220      gsl_odeiv2_driver_free (state);
221      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
222    }
223
224  /* Distribute pointer to driver object */
225
226  gsl_odeiv2_step_set_driver (state->s, state);
227  gsl_odeiv2_evolve_set_driver (state->e, state);
228  gsl_odeiv2_control_set_driver (state->c, state);
229
230  return state;
231}
232
233gsl_odeiv2_driver *
234gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * sys,
235                                      const gsl_odeiv2_step_type * T,
236                                      const double hstart,
237                                      const double epsabs,
238                                      const double epsrel, const double a_y,
239                                      const double a_dydt)
240{
241  /* Initializes an ODE driver system with control object of type
242     standard_new. 
243   */
244
245  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
246
247  if (state == NULL)
248    {
249      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
250    }
251
252  if (epsabs >= 0.0 && epsrel >= 0.0)
253    {
254      state->c =
255        gsl_odeiv2_control_standard_new (epsabs, epsrel, a_y, a_dydt);
256
257      if (state->c == NULL)
258        {
259          gsl_odeiv2_driver_free (state);
260          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
261        }
262    }
263  else
264    {
265      gsl_odeiv2_driver_free (state);
266      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
267    }
268
269  /* Distribute pointer to driver object */
270
271  gsl_odeiv2_step_set_driver (state->s, state);
272  gsl_odeiv2_evolve_set_driver (state->e, state);
273  gsl_odeiv2_control_set_driver (state->c, state);
274
275  return state;
276}
277
278gsl_odeiv2_driver *
279gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * sys,
280                                    const gsl_odeiv2_step_type * T,
281                                    const double hstart,
282                                    const double epsabs, const double epsrel,
283                                    const double a_y, const double a_dydt,
284                                    const double scale_abs[])
285{
286  /* Initializes an ODE driver system with control object of type
287     scaled_new. 
288   */
289
290  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);
291
292  if (state == NULL)
293    {
294      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
295    }
296
297  if (epsabs >= 0.0 && epsrel >= 0.0)
298    {
299      state->c = gsl_odeiv2_control_scaled_new (epsabs, epsrel, a_y, a_dydt,
300                                                scale_abs, sys->dimension);
301
302      if (state->c == NULL)
303        {
304          gsl_odeiv2_driver_free (state);
305          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
306        }
307    }
308  else
309    {
310      gsl_odeiv2_driver_free (state);
311      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
312    }
313
314  /* Distribute pointer to driver object */
315
316  gsl_odeiv2_step_set_driver (state->s, state);
317  gsl_odeiv2_evolve_set_driver (state->e, state);
318  gsl_odeiv2_control_set_driver (state->c, state);
319
320  return state;
321}
322
323int
324gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t,
325                         const double t1, double y[])
326{
327  /* Main driver function that evolves the system from t to t1. In
328     beginning vector y contains the values of dependent variables at
329     t. This function returns values at t=t1 in y. In case of
330     unrecoverable error, y and t contains the values after the last
331     successful step.
332   */
333
334  int sign = 0;
335  d->n = 0;
336
337  /* Determine integration direction sign */
338
339  if (d->h > 0.0)
340    {
341      sign = 1;
342    }
343  else
344    {
345      sign = -1;
346    }
347
348  /* Check that t, t1 and step direction are sensible */
349
350  if (sign * (t1 - *t) < 0.0)
351    {
352      GSL_ERROR_NULL
353        ("integration limits and/or step direction not consistent",
354         GSL_EINVAL);
355    }
356
357  /* Evolution loop */
358
359  while (sign * (t1 - *t) > 0.0)
360    {
361      int s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, d->sys,
362                                       t, t1, &(d->h), y);
363
364      if (s != GSL_SUCCESS)
365        {
366          return s;
367        }
368
369      /* Check for maximum allowed steps */
370
371      if ((d->nmax > 0) && (d->n > d->nmax))
372        {
373          return GSL_EMAXITER;
374        }
375
376      /* Set step size if maximum size is exceeded */
377
378      if (fabs (d->h) > d->hmax)
379        {
380          d->h = sign * d->hmax;
381        }
382
383      /* Check for too small step size */
384
385      if (fabs (d->h) < d->hmin)
386        {
387          return GSL_ENOPROG;
388        }
389
390      d->n++;
391    }
392
393  return GSL_SUCCESS;
394}
395
396int
397gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t,
398                                    const double h, const unsigned long int n,
399                                    double y[])
400{
401  /* Alternative driver function that evolves the system from t using
402   * n steps of size h. In the beginning vector y contains the values
403   * of dependent variables at t. This function returns values at t =
404   * t + n * h in y. In case of an unrecoverable error, y and t
405   * contains the values after the last successful step.
406   */
407
408  unsigned long int i;
409  d->n = 0;
410
411  /* Evolution loop */
412
413  for (i = 0; i < n; i++)
414    {
415      int s = gsl_odeiv2_evolve_apply_fixed_step (d->e, d->c, d->s, d->sys,
416                                                  t, h, y);
417
418      if (s != GSL_SUCCESS)
419        {
420          return s;
421        }
422
423      d->n++;
424    }
425
426  return GSL_SUCCESS;
427}
428
429int
430gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d)
431{
432  /* Reset the driver. Resets evolve and step objects. */
433
434  {
435    int s = gsl_odeiv2_evolve_reset (d->e);
436
437    if (s != GSL_SUCCESS)
438      {
439        return s;
440      }
441  }
442
443  {
444    int s = gsl_odeiv2_step_reset (d->s);
445
446    if (s != GSL_SUCCESS)
447      {
448        return s;
449      }
450  }
451
452  return GSL_SUCCESS;
453}
454
455int
456gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart)
457{
458  /* Resets current driver and sets initial step size to hstart */
459
460  gsl_odeiv2_driver_reset (d);
461
462  if ((d->hmin > fabs (hstart)) || (fabs (hstart) > d->hmax))
463    {
464      GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
465    }
466
467  if (hstart > 0.0 || hstart < 0.0)
468    {
469      d->h = hstart;
470    }
471  else
472    {
473      GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL);
474    }
475
476  return GSL_SUCCESS;
477}
478
479void
480gsl_odeiv2_driver_free (gsl_odeiv2_driver * state)
481{
482  if (state->c != NULL)
483    {
484      gsl_odeiv2_control_free (state->c);
485    }
486
487  gsl_odeiv2_evolve_free (state->e);
488  gsl_odeiv2_step_free (state->s);
489  free (state);
490}
491
Events list
Event 1
Event 2
Event 3
Event 4
Event 5
Event 6
Event 12
Event 13
Event 26
Event 27
Event 28
Event 29