Annotated Source Code

1/* ode-initval/evolve.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/* Author:  G. Jungman
21 */
22#include <config.h>
23#include <string.h>
24#include <stdlib.h>
25#include <gsl/gsl_math.h>
26#include <gsl/gsl_errno.h>
27#include <gsl/gsl_odeiv2.h>
28
29#include "odeiv_util.h"
30
31gsl_odeiv2_evolve *
32gsl_odeiv2_evolve_alloc (size_t dim)
14Enter function.    gsl_odeiv2_evolve_alloc(dim)
33{
34  gsl_odeiv2_evolve *e =
35    (gsl_odeiv2_evolve *) malloc (sizeof (gsl_odeiv2_evolve));
15Call a function.    malloc(sizeof(gsl_odeiv2_evolve))
36
37  if (e == 0)
16Take the false branch.    e == 0
38    {
39      GSL_ERROR_NULL ("failed to allocate space for evolve struct",
40                      GSL_ENOMEM);
41    }
42
43  e->y0 = (double *) malloc (dim * sizeof (double));
17Call a function.    malloc(dim * sizeof(double))
44
45  if (e->y0 == 0)
18Take the false branch.    e->y0 == 0
46    {
47      free (e);
48      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
49    }
50
51  e->yerr = (double *) malloc (dim * sizeof (double));
19Call a function.    malloc(dim * sizeof(double))
52
53  if (e->yerr == 0)
20Take the false branch.    e->yerr == 0
54    {
55      free (e->y0);
56      free (e);
57      GSL_ERROR_NULL ("failed to allocate space for yerr", GSL_ENOMEM);
58    }
59
60  e->dydt_in = (double *) malloc (dim * sizeof (double));
21Call a function.    malloc(dim * sizeof(double))
33Memory allocated in heap space is not freed.    malloc(dim * sizeof(double))
61
62  if (e->dydt_in == 0)
22Take the false branch.    e->dydt_in == 0
63    {
64      free (e->yerr);
65      free (e->y0);
66      free (e);
67      GSL_ERROR_NULL ("failed to allocate space for dydt_in", GSL_ENOMEM);
68    }
69
70  e->dydt_out = (double *) malloc (dim * sizeof (double));
23Call a function.    malloc(dim * sizeof(double))
71
72  if (e->dydt_out == 0)
24Take the false branch.    e->dydt_out == 0
73    {
74      free (e->dydt_in);
75      free (e->yerr);
76      free (e->y0);
77      free (e);
78      GSL_ERROR_NULL ("failed to allocate space for dydt_out", GSL_ENOMEM);
79    }
80
81  e->dimension = dim;
82  e->count = 0;
83  e->failed_steps = 0;
84  e->last_step = 0.0;
85  e->driver = NULL;
86
87  return e;
88}
25Exit function.    gsl_odeiv2_evolve_alloc(dim)
89
90int
91gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * e)
92{
93  e->count = 0;
94  e->failed_steps = 0;
95  e->last_step = 0.0;
96  return GSL_SUCCESS;
97}
98
99void
100gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * e)
101{
102  RETURN_IF_NULL (e);
103  free (e->dydt_out);
104  free (e->dydt_in);
105  free (e->yerr);
106  free (e->y0);
107  free (e);
108}
109
110/* Evolution framework method.
111 *
112 * Uses an adaptive step control object
113 */
114int
115gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * e,
116                         gsl_odeiv2_control * con,
117                         gsl_odeiv2_step * step,
118                         const gsl_odeiv2_system * dydt,
119                         double *t, double t1, double *h, double y[])
120{
121  const double t0 = *t;
122  double h0 = *h;
123  int step_status;
124  int final_step = 0;
125  double dt = t1 - t0;          /* remaining time, possibly less than h */
126
127  if (e->dimension != step->dimension)
128    {
129      GSL_ERROR ("step dimension must match evolution size", GSL_EINVAL);
130    }
131
132  if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0))
133    {
134      GSL_ERROR ("step direction must match interval direction", GSL_EINVAL);
135    }
136
137  /* Save y in case of failure in a step */
138
139  DBL_MEMCPY (e->y0, y, e->dimension);
140
141  /* Calculate initial dydt once or reuse previous value if the method
142     can benefit. */
143  
144  if (step->type->can_use_dydt_in)
145    {
146      if (e->count == 0)
147        {
148          int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in);
149
150          if (status)
151            {
152              return status;
153            }
154        }
155      else
156        {
157          DBL_MEMCPY (e->dydt_in, e->dydt_out, e->dimension);
158        }
159    }
160
161try_step:
162
163  if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt))
164    {
165      h0 = dt;
166      final_step = 1;
167    }
168  else
169    {
170      final_step = 0;
171    }
172
173  if (step->type->can_use_dydt_in)
174    {
175      step_status =
176        gsl_odeiv2_step_apply (step, t0, h0, y, e->yerr, e->dydt_in,
177                               e->dydt_out, dydt);
178    }
179  else
180    {
181      step_status =
182        gsl_odeiv2_step_apply (step, t0, h0, y, e->yerr, NULL, e->dydt_out,
183                               dydt);
184    }
185
186  /* Return if stepper indicates a pointer or user function failure */
187
188  if (step_status == GSL_EFAULT || step_status == GSL_EBADFUNC)
189    {
190      return step_status;
191    }
192
193  /* Check for stepper internal failure */
194
195  if (step_status != GSL_SUCCESS)
196    {
197      /* Stepper was unable to calculate step. Try decreasing step size. */
198
199      double h_old = h0;
200
201      h0 *= 0.5;
202
203      /* Check that an actual decrease in h0 occured and the
204         suggested h0 will change the time by at least 1 ulp */
205
206      {
207        double t_curr = GSL_COERCE_DBL (*t);
208        double t_next = GSL_COERCE_DBL ((*t) + h0);
209
210        if (fabs (h0) < fabs (h_old) && t_next != t_curr)
211          {
212            
213            /* Step was decreased. Undo step, and try again with new h0. */
214            
215            DBL_MEMCPY (y, e->y0, dydt->dimension);
216            e->failed_steps++;
217            goto try_step;
218          }
219        else
220          {
221            *h = h0;              /* notify user of step-size which caused the failure */
222            *t = t0;              /* restore original t value */
223            return step_status;
224          }
225      }
226    }
227
228  e->count++;
229  e->last_step = h0;
230
231  if (final_step)
232    {
233      *t = t1;
234    }
235  else
236    {
237      *t = t0 + h0;
238    }
239
240  if (con != NULL)
241    {
242      /* Check error and attempt to adjust the step. */
243
244      double h_old = h0;
245
246      const int hadjust_status
247        =
248        gsl_odeiv2_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);
249
250      if (hadjust_status == GSL_ODEIV_HADJ_DEC)
251        {
252          /* Check that the reported status is correct (i.e. an actual
253             decrease in h0 occured) and the suggested h0 will change
254             the time by at least 1 ulp */
255
256          double t_curr = GSL_COERCE_DBL (*t);
257          double t_next = GSL_COERCE_DBL ((*t) + h0);
258
259          if (fabs (h0) < fabs (h_old) && t_next != t_curr)
260            {
261              /* Step was decreased. Undo step, and try again with new h0. */
262
263              DBL_MEMCPY (y, e->y0, dydt->dimension);
264              e->failed_steps++;
265              goto try_step;
266            }
267          else
268            {
269              /* Can not obtain required error tolerance, and can not
270                 decrease step-size any further, so give up and return
271                 GSL_FAILURE.
272               */
273
274              *h = h0;          /* notify user of step-size which caused the failure */
275              return GSL_FAILURE;
276            }
277        }
278    }
279
280  /* Suggest step size for next time-step. Change of step size is not
281     suggested in the final step, because that step can be very
282     small compared to previous step, to reach t1. 
283   */
284
285  if (final_step == 0)
286    {
287      *h = h0;
288    }
289
290  return step_status;
291}
292
293/* Evolves the system using the user specified constant step size h.
294 */
295
296int
297gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * e,
298                                    gsl_odeiv2_control * con,
299                                    gsl_odeiv2_step * step,
300                                    const gsl_odeiv2_system * dydt,
301                                    double *t, const double h, double y[])
302{
303  const double t0 = *t;
304  int step_status;
305
306  if (e->dimension != step->dimension)
307    {
308      GSL_ERROR ("step dimension must match evolution size", GSL_EINVAL);
309    }
310
311  /* Save y in case of failure in a step */
312
313  DBL_MEMCPY (e->y0, y, e->dimension);
314
315  /* Calculate initial dydt once if the method can benefit. */
316
317  if (step->type->can_use_dydt_in)
318    {
319      int status = GSL_ODEIV_FN_EVAL (dydt, t0, y, e->dydt_in);
320
321      if (status)
322        {
323          return status;
324        }
325    }
326
327  if (step->type->can_use_dydt_in)
328    {
329      step_status =
330        gsl_odeiv2_step_apply (step, t0, h, y, e->yerr, e->dydt_in,
331                               e->dydt_out, dydt);
332    }
333  else
334    {
335      step_status =
336        gsl_odeiv2_step_apply (step, t0, h, y, e->yerr, NULL, e->dydt_out,
337                               dydt);
338    }
339
340  /* Return the stepper return value in case of an error */
341
342  if (step_status != GSL_SUCCESS)
343    {
344      return step_status;
345    }
346
347  if (con != NULL)
348    {
349      /* Calculate error level. Fail if error level exceeds desired
350         error level. */
351
352      double htemp = h;
353
354      const int hadjust_status
355        = gsl_odeiv2_control_hadjust (con, step, y, e->yerr,
356                                      e->dydt_out, &htemp);
357
358      if (hadjust_status == GSL_ODEIV_HADJ_DEC)
359        {
360          DBL_MEMCPY (y, e->y0, dydt->dimension);
361          e->failed_steps++;
362          return GSL_FAILURE;
363        }
364    }
365
366  /* Step is accepted, update status */
367
368  e->count++;
369  e->last_step = h;
370  *t = t0 + h;
371
372  return GSL_SUCCESS;
373}
374
375int
376gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * e,
377                              const gsl_odeiv2_driver * d)
378{
379  if (d != NULL)
380    {
381      e->driver = d;
382    }
383  else
384    {
385      GSL_ERROR_NULL ("driver pointer is null", GSL_EFAULT);
386    }
387
388  return GSL_SUCCESS;
389}
390
Events list
Event 14
Event 15
Event 16
Event 17
Event 18
Event 19
Event 20
Event 21
Event 33
Event 22
Event 23
Event 24
Event 25