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 | |
31 | gsl_odeiv2_evolve * |
32 | gsl_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)) |
| 33Memory allocated in heap space is not freed.    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)) |
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 | |
90 | int |
91 | gsl_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 | |
99 | void |
100 | gsl_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 | */ |
114 | int |
115 | gsl_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 | |
161 | try_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 | |
296 | int |
297 | gsl_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 | |
375 | int |
376 | gsl_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 | |