Bulirsch-Stoer method implementation for solving ordinary differential equations using polynomial extrapolation. More...
#include "mdb.h"Go to the source code of this file.
Macros | |
| #define | DEBUG 0 |
| #define | IMAX 11 |
| #define | NUSE 7 |
| #define | TINY 1.0e-30 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
| #define | TINY 1.0e-30 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
Functions | |
| void | new_scale_factors_dp (double *yscale, double *y0, double *dydx0, double h_start, double *tiny, long *accmode, double *accuracy, long n_eq) |
| void | initial_scale_factors_dp (double *yscale, double *y0, double *dydx0, double h_start, double *tiny, long *accmode, double *accuracy, double *accur, double x0, double xf, long n_eq) |
| void | bs_qctune (double newStepIncreaseFactor, double newStepDecreaseFactor) |
| long | bs_step (double *yFinal, double *x, double *yInitial, double *dydxInitial, double step, double *stepUsed, double *stepRecommended, double *yScale, long equations, void(*derivs)(double *dydx, double *y, double x), long *misses) |
| long | bs_odeint (double *y0, void(*derivs)(double *dydx, double *y, double x), long n_eq, double *accuracy, long *accmode, double *tiny, long *misses, double *x0, double xf, double x_accuracy, double h_start, double h_max, double *h_rec, double(*exit_func)(double *dydx, double *y, double x), double exit_accuracy, long n_to_skip, void(*store_data)(double *dydx, double *y, double x, double exf)) |
| Integrates a system of ordinary differential equations using the Bulirsch-Stoer method. | |
| long | bs_odeint1 (double *y0, void(*derivs)(), long n_eq, double *accuracy, long *accmode, double *tiny, long *misses, double *x0, double xf, double x_accuracy, double h_start, double h_max, double *h_rec) |
| Integrates a system of ordinary differential equations to a specified upper limit without exit condition checks or intermediate data storage. | |
| long | bs_odeint2 (double *y0, void(*derivs)(double *dydx, double *y, double x), long n_eq, double *accuracy, long *accmode, double *tiny, long *misses, double *x0, double xf, double x_accuracy, double h_start, double h_max, double *h_rec, double exit_value, long i_exit_value, double exit_accuracy, long n_to_skip) |
| Integrates a system of ordinary differential equations until a specified component reaches a target value or the upper limit is met. | |
| long | bs_odeint3 (double *yif, void(*derivs)(double *dydx, double *y, double x), long n_eq, double *accuracy, long *accmode, double *tiny, long *misses, double *x0, double xf, double x_accuracy, double h_start, double h_max, double *h_rec, double(*exit_func)(double *dydx, double *y, double x), double exit_accuracy) |
| Integrates a system of ordinary differential equations using the Bulirsch-Stoer method with optimized internal state management. | |
| long | bs_odeint4 (double *y0, void(*derivs)(double *dydx, double *y, double x), long n_eq, double *accuracy, long *accmode, double *tiny, long *misses, double *x0, double xf, double x_accuracy, double h_start, double h_max, double *h_rec, double exit_value, long i_exit_value, double exit_accuracy, long n_to_skip, void(*store_data)(double *dydx, double *y, double x, double exf)) |
| Integrates a system of ordinary differential equations until a specified component reaches a target value or the upper limit is met, with intermediate data storage. | |
Variables | |
| static double | stepIncreaseFactor = 0.50 |
| static double | stepDecreaseFactor = 0.95 |
Detailed Description
Bulirsch-Stoer method implementation for solving ordinary differential equations using polynomial extrapolation.
This file contains the implementation of the Bulirsch-Stoer method for integrating ordinary differential equations (ODEs). It includes functions for performing integration steps, handling scale factors, and managing accuracy controls.
- Copyright
- (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
- (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
- License
- This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.
Definition in file bsODEp.c.
Macro Definition Documentation
◆ DEBUG
◆ IMAX
◆ NUSE
◆ TINY [1/2]
◆ TINY [2/2]
Function Documentation
◆ bs_odeint()
| long bs_odeint | ( | double * | y0, |
| void(* | derivs )(double *dydx, double *y, double x), | ||
| long | n_eq, | ||
| double * | accuracy, | ||
| long * | accmode, | ||
| double * | tiny, | ||
| long * | misses, | ||
| double * | x0, | ||
| double | xf, | ||
| double | x_accuracy, | ||
| double | h_start, | ||
| double | h_max, | ||
| double * | h_rec, | ||
| double(* | exit_func )(double *dydx, double *y, double x), | ||
| double | exit_accuracy, | ||
| long | n_to_skip, | ||
| void(* | store_data )(double *dydx, double *y, double x, double exf) ) |
Integrates a system of ordinary differential equations using the Bulirsch-Stoer method.
This function integrates a set of ODEs from an initial value until the upper limit of the independent variable is reached or a user-supplied exit condition function evaluates to zero.
- Parameters
-
y0 Pointer to the array of initial values of the dependent variables. Upon successful completion, it contains the final values. derivs Function pointer to compute the derivatives. It calculates dy/dx given the current state. n_eq Number of equations or dependent variables in the system. accuracy Pointer to the array specifying the desired accuracy for each dependent variable. accmode Pointer to the array specifying the accuracy control mode for each dependent variable. Modes: - 0: Fractional accuracy per step for each variable.
- 1: Fractional accuracy globally.
- 2: Absolute accuracy per step for each variable.
- 3: Absolute accuracy globally.
tiny Pointer to the array of small values representing the lower limits of significance for each dependent variable. misses Pointer to the array that counts the number of times each variable caused a step size reset due to exceeding error tolerance. x0 Pointer to the initial value of the independent variable. It is updated to the final value after integration. xf Upper limit of the independent variable to integrate up to. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size for the integration. h_max Maximum allowed step size for the integration. h_rec Pointer to the variable where the recommended step size for continuation will be stored. exit_func Function pointer to the exit condition function. It returns a value that, when zero, signals the integration to stop. exit_accuracy Desired accuracy for the exit condition function to evaluate to zero. n_to_skip Number of zeros of the exit function to skip before returning. store_data Function pointer to store intermediate integration points. It is called with the current derivatives, dependent variables, independent variable, and exit function value.
- Returns
- Returns a non-negative value on failure (with specific codes) or a positive value on success.
Definition at line 179 of file bsODEp.c.
203 {
204 double *y_return, *accur;
205 double *dydx0, *y1, *dydx1, *dydx2, *y2;
206 double ex0, ex1, ex2, x1, x2, *yscale;
207 double h_used, h_next, xdiff;
208 long i, n_step_ups = 0, is_zero;
209#define MAX_N_STEP_UPS 10
210
211 if (*x0 > xf)
212 return (DIFFEQ_XI_GT_XF);
213 if (FABS(*x0 - xf) < x_accuracy)
214 return (DIFFEQ_SOLVED_ALREADY);
215
216 /* Meaning of accmode:
217 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
218 * each step for ith variable. tiny[i] is lower limit
219 * of significance for the ith variable.
220 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
221 * satisfied globally, not locally.
222 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
223 * step for the ith variable. tiny[i] is ignored.
224 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
225 * be satisfied globally, not locally.
226 */
227 for (i = 0; i < n_eq; i++) {
228 if (accmode[i] < 0 || accmode[i] > 3)
230 if (accmode[i] < 2 && tiny[i] < TINY)
231 tiny[i] = TINY;
232 misses[i] = 0;
233 }
234
235 y_return = y0;
242
243 /* calculate derivatives and exit function at the initial point */
244 (*derivs)(dydx0, y0, *x0);
245
246 /* set the scales for evaluating accuracy. yscale[i] is the
247 * absolute level of accuracy required of the next integration step
248 */
250 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
251 accuracy, accur, *x0, xf, n_eq);
252
253 ex0 = exit_func ? (*exit_func)(dydx0, y0, *x0) : 0;
254 if (store_data)
255 (*store_data)(dydx0, y0, *x0, ex0);
256 is_zero = 0;
257
258 do {
259 /* check for zero of exit function */
260 if (exit_func && FABS(ex0) < exit_accuracy) {
261 if (!is_zero) {
262 if (n_to_skip == 0) {
263 if (store_data)
264 (*store_data)(dydx0, y0, *x0, ex0);
265 for (i = 0; i < n_eq; i++)
266 y_return[i] = y0[i];
267 *h_rec = h_start;
268 tfree(dydx0);
269 tfree(dydx1);
270 tfree(dydx2);
271 tfree(yscale);
272 tfree(accur);
273 if (y0 != y_return)
274 tfree(y0);
275 if (y1 != y_return)
276 tfree(y1);
277 if (y2 != y_return)
278 tfree(y2);
279 return (DIFFEQ_ZERO_FOUND);
280 } else {
281 is_zero = 1;
282 --n_to_skip;
283 }
284 }
285 } else
286 is_zero = 0;
287 /* adjust step size to stay within interval */
288 if ((xdiff = xf - *x0) < h_start)
289 h_start = xdiff;
290 /* take a step */
291 x1 = *x0;
292 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
293 yscale, n_eq, derivs, misses)) {
294 if (n_step_ups++ > MAX_N_STEP_UPS)
296 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
297 continue;
298 }
299 /* calculate derivatives and exit function at new point */
300 (*derivs)(dydx1, y1, x1);
301 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
302 if (store_data)
303 (*store_data)(dydx1, y1, x1, ex1);
304 /* check for change in sign of exit function */
305 if (exit_func && SIGN(ex0) != SIGN(ex1) && !is_zero) {
306 if (n_to_skip == 0)
307 break;
308 else {
309 --n_to_skip;
310 is_zero = 1;
311 }
312 }
313 /* check for end of interval */
314 if (FABS(xdiff = xf - x1) < x_accuracy) {
315 /* end of the interval */
316 if (store_data) {
317 (*derivs)(dydx1, y1, x1);
318 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
319 (*store_data)(dydx1, y1, x1, ex1);
320 }
321 for (i = 0; i < n_eq; i++)
322 y_return[i] = y1[i];
323 *x0 = x1;
324 *h_rec = h_start;
325 tfree(dydx0);
326 tfree(dydx1);
327 tfree(dydx2);
328 tfree(yscale);
329 tfree(accur);
330 if (y0 != y_return)
331 tfree(y0);
332 if (y1 != y_return)
333 tfree(y1);
334 if (y2 != y_return)
335 tfree(y2);
336 return (DIFFEQ_END_OF_INTERVAL);
337 }
338 /* copy the new solution into the old variables */
339 SWAP_PTR(dydx0, dydx1);
340 SWAP_PTR(y0, y1);
341 ex0 = ex1;
342 *x0 = x1;
343 /* adjust the step size as recommended by bs_step() */
344 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
345 /* calculate new scale factors */
346 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
347 accur, n_eq);
348 } while (1);
349 *h_rec = h_start;
350
351 if (!exit_func) {
352 printf("failure in bs_odeint(): solution stepped outside interval\n");
353 tfree(dydx0);
354 tfree(dydx1);
355 tfree(dydx2);
356 tfree(yscale);
357 tfree(accur);
358 if (y0 != y_return)
359 tfree(y0);
360 if (y1 != y_return)
361 tfree(y1);
362 if (y2 != y_return)
363 tfree(y2);
364 return (DIFFEQ_OUTSIDE_INTERVAL);
365 }
366
367 if (FABS(ex1) < exit_accuracy) {
368 for (i = 0; i < n_eq; i++)
369 y_return[i] = y1[i];
370 *x0 = x1;
371 tfree(dydx0);
372 tfree(dydx1);
373 tfree(dydx2);
374 tfree(yscale);
375 tfree(accur);
376 if (y0 != y_return)
377 tfree(y0);
378 if (y1 != y_return)
379 tfree(y1);
380 if (y2 != y_return)
381 tfree(y2);
382 return (DIFFEQ_ZERO_FOUND);
383 }
384
385 /* The root has been bracketed. */
386 do {
387 /* try to take a step to the position where the zero is expected */
388 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
389 x2 = *x0;
390 /* calculate new scale factors */
391 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
392 accur, n_eq);
393 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
394 yscale, n_eq, derivs, misses))
396 /* check the exit function at the new position */
397 (*derivs)(dydx2, y2, x2);
398 ex2 = (*exit_func)(dydx2, y2, x2);
399 if (FABS(ex2) < exit_accuracy) {
400 for (i = 0; i < n_eq; i++)
401 y_return[i] = y2[i];
402 *x0 = x2;
403 tfree(dydx0);
404 tfree(dydx1);
405 tfree(dydx2);
406 tfree(yscale);
407 tfree(accur);
408 if (y0 != y_return)
409 tfree(y0);
410 if (y1 != y_return)
411 tfree(y1);
412 if (y2 != y_return)
413 tfree(y2);
414 return (DIFFEQ_ZERO_FOUND);
415 }
416 /* rebracket the root */
417 if (SIGN(ex1) == SIGN(ex2)) {
418 SWAP_PTR(y1, y2);
419 SWAP_PTR(dydx1, dydx2);
420 x1 = x2;
421 ex1 = ex2;
422 } else {
423 SWAP_PTR(y0, y2);
424 SWAP_PTR(dydx0, dydx2);
425 *x0 = x2;
426 ex0 = ex2;
427 }
428 } while (1);
429}
int tfree(void *ptr)
Frees a memory block and records the deallocation if tracking is enabled.
Definition array.c:230
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
◆ bs_odeint1()
| long bs_odeint1 | ( | double * | y0, |
| void(* | derivs )(), | ||
| long | n_eq, | ||
| double * | accuracy, | ||
| long * | accmode, | ||
| double * | tiny, | ||
| long * | misses, | ||
| double * | x0, | ||
| double | xf, | ||
| double | x_accuracy, | ||
| double | h_start, | ||
| double | h_max, | ||
| double * | h_rec ) |
Integrates a system of ordinary differential equations to a specified upper limit without exit condition checks or intermediate data storage.
This function is a streamlined version of bs_odeint() that integrates the ODE system until the upper limit of the independent variable is reached. It does not evaluate a user-supplied exit condition function or store intermediate integration points, resulting in improved performance.
- Parameters
-
y0 Pointer to the array of initial values of the dependent variables. Upon successful completion, it contains the final values. derivs Function pointer to compute the derivatives. It calculates dy/dx given the current state. n_eq Number of equations or dependent variables in the system. accuracy Pointer to the array specifying the desired accuracy for each dependent variable. accmode Pointer to the array specifying the accuracy control mode for each dependent variable. Modes: - 0: Fractional accuracy per step for each variable.
- 1: Fractional accuracy globally.
- 2: Absolute accuracy per step for each variable.
- 3: Absolute accuracy globally.
tiny Pointer to the array of small values representing the lower limits of significance for each dependent variable. misses Pointer to the array that counts the number of times each variable caused a step size reset due to exceeding error tolerance. x0 Pointer to the initial value of the independent variable. It is updated to the final value after integration. xf Upper limit of the independent variable to integrate up to. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size for the integration. h_max Maximum allowed step size for the integration. h_rec Pointer to the variable where the recommended step size for continuation will be stored.
- Returns
- Returns 1 on successful integration, or a non-negative value on failure.
Definition at line 456 of file bsODEp.c.
473 {
474 double *y_return;
475 double *dydx0, *y1, *dydx1, *dydx2, *y2;
476 double x1, *yscale, *accur;
477 double h_used, h_next, xdiff;
478 long i, n_step_ups = 0;
479#define MAX_N_STEP_UPS 10
480
481 if (*x0 > xf)
482 return (DIFFEQ_XI_GT_XF);
483 if (fabs(*x0 - xf) < x_accuracy)
484 return (DIFFEQ_SOLVED_ALREADY);
485
486 /* Meaning of accmode:
487 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
488 * each step for ith variable. tiny[i] is lower limit
489 * of significance for the ith variable.
490 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
491 * satisfied globally, not locally.
492 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
493 * step for the ith variable. tiny[i] is ignored.
494 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
495 * be satisfied globally, not locally.
496 */
497 for (i = 0; i < n_eq; i++) {
498 if (accmode[i] < 0 || accmode[i] > 3)
500 if (accmode[i] < 2 && tiny[i] < TINY)
501 tiny[i] = TINY;
502 misses[i] = 0;
503 }
504
505 y_return = y0;
512
513 /* calculate derivatives at the initial point */
514 (*derivs)(dydx0, y0, *x0);
515
516 /* set the scales for evaluating accuracy. yscale[i] is the
517 * absolute level of accuracy required of the next integration step
518 */
520 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
521 accuracy, accur, *x0, xf, n_eq);
522
523 do {
524 /* adjust step size to stay within interval */
525 if ((xdiff = xf - *x0) < h_start)
526 h_start = xdiff;
527 /* take a step */
528 x1 = *x0;
529 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
530 yscale, n_eq, derivs, misses)) {
531 if (n_step_ups++ > MAX_N_STEP_UPS)
533 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
534 continue;
535 }
536 /* check for end of interval */
537 if (fabs(xdiff = xf - x1) < x_accuracy) {
538 /* end of the interval */
539 for (i = 0; i < n_eq; i++)
540 y_return[i] = y1[i];
541 *x0 = x1;
542 *h_rec = h_start;
543 tfree(dydx0);
544 tfree(dydx1);
545 tfree(dydx2);
546 tfree(yscale);
547 tfree(accur);
548 if (y0 != y_return)
549 tfree(y0);
550 if (y1 != y_return)
551 tfree(y1);
552 if (y2 != y_return)
553 tfree(y2);
554 return (DIFFEQ_END_OF_INTERVAL);
555 }
556 /* calculate derivatives at new point */
557 (*derivs)(dydx1, y1, x1);
558 /* copy the new solution into the old variables */
559 SWAP_PTR(dydx0, dydx1);
560 SWAP_PTR(y0, y1);
561 *x0 = x1;
562 /* adjust the step size as recommended by bs_step() */
563 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
564 /* calculate new scale factors */
565 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
566 accur, n_eq);
567 } while (1);
568}
◆ bs_odeint2()
| long bs_odeint2 | ( | double * | y0, |
| void(* | derivs )(double *dydx, double *y, double x), | ||
| long | n_eq, | ||
| double * | accuracy, | ||
| long * | accmode, | ||
| double * | tiny, | ||
| long * | misses, | ||
| double * | x0, | ||
| double | xf, | ||
| double | x_accuracy, | ||
| double | h_start, | ||
| double | h_max, | ||
| double * | h_rec, | ||
| double | exit_value, | ||
| long | i_exit_value, | ||
| double | exit_accuracy, | ||
| long | n_to_skip ) |
Integrates a system of ordinary differential equations until a specified component reaches a target value or the upper limit is met.
This function integrates the ODE system until either the independent variable reaches the specified upper limit or a particular dependent variable reaches a target value within a specified accuracy. It does not evaluate a general exit condition function or store intermediate data, offering improved performance compared to bs_odeint().
- Parameters
-
y0 Pointer to the array of initial values of the dependent variables. Upon successful completion, it contains the final values. derivs Function pointer to compute the derivatives. It calculates dy/dx given the current state. n_eq Number of equations or dependent variables in the system. accuracy Pointer to the array specifying the desired accuracy for each dependent variable. accmode Pointer to the array specifying the accuracy control mode for each dependent variable. Modes: - 0: Fractional accuracy per step for each variable.
- 1: Fractional accuracy globally.
- 2: Absolute accuracy per step for each variable.
- 3: Absolute accuracy globally.
tiny Pointer to the array of small values representing the lower limits of significance for each dependent variable. misses Pointer to the array that counts the number of times each variable caused a step size reset due to exceeding error tolerance. x0 Pointer to the initial value of the independent variable. It is updated to the final value after integration. xf Upper limit of the independent variable to integrate up to. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size for the integration. h_max Maximum allowed step size for the integration. h_rec Pointer to the variable where the recommended step size for continuation will be stored. exit_value Target value that the specified component of the solution should reach. i_exit_value Index of the dependent variable component that is being monitored for reaching the target value. exit_accuracy Desired accuracy for the target value condition. n_to_skip Number of times the target condition can be met before integration stops.
- Returns
- Returns 1 on successful integration, or a non-negative value on failure.
Definition at line 599 of file bsODEp.c.
622 {
623 double *y_return, *accur;
624 double *dydx0, *y1, *dydx1, *dydx2, *y2;
625 double ex0, ex1, ex2, x1, x2, *yscale;
626 double h_used, h_next, xdiff;
627 long i, n_step_ups = 0, is_zero;
628#define MAX_N_STEP_UPS 10
629
630 if (*x0 > xf)
631 return (DIFFEQ_XI_GT_XF);
632 if (fabs(*x0 - xf) < x_accuracy)
633 return (DIFFEQ_SOLVED_ALREADY);
634 if (i_exit_value < 0 || i_exit_value >= n_eq)
636
637 /* Meaning of accmode:
638 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
639 * each step for ith variable. tiny[i] is lower limit
640 * of significance for the ith variable.
641 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
642 * satisfied globally, not locally.
643 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
644 * step for the ith variable. tiny[i] is ignored.
645 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
646 * be satisfied globally, not locally.
647 */
648 for (i = 0; i < n_eq; i++) {
649 if (accmode[i] < 0 || accmode[i] > 3)
651 if (accmode[i] < 2 && tiny[i] < TINY)
652 tiny[i] = TINY;
653 misses[i] = 0;
654 }
655
656 y_return = y0;
663
664 /* calculate derivatives and exit function at the initial point */
665 (*derivs)(dydx0, y0, *x0);
666
667 /* set the scales for evaluating accuracy. yscale[i] is the
668 * absolute level of accuracy required of the next integration step
669 */
671 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
672 accuracy, accur, *x0, xf, n_eq);
673
674 ex0 = exit_value - y0[i_exit_value];
675 is_zero = 0;
676 do {
677 /* check for zero of exit function */
678 if (fabs(ex0) < exit_accuracy) {
679 if (!is_zero) {
680 if (n_to_skip == 0) {
681 for (i = 0; i < n_eq; i++)
682 y_return[i] = y0[i];
683 *h_rec = h_start;
684 tfree(dydx0);
685 tfree(dydx1);
686 tfree(dydx2);
687 tfree(yscale);
688 tfree(accur);
689 if (y0 != y_return)
690 tfree(y0);
691 if (y1 != y_return)
692 tfree(y1);
693 if (y2 != y_return)
694 tfree(y2);
695 return (DIFFEQ_ZERO_FOUND);
696 } else {
697 is_zero = 1;
698 --n_to_skip;
699 }
700 }
701 } else
702 is_zero = 0;
703 /* adjust step size to stay within interval */
704 if ((xdiff = xf - *x0) < h_start)
705 h_start = xdiff;
706 /* take a step */
707 x1 = *x0;
708 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
709 yscale, n_eq, derivs, misses)) {
710 if (n_step_ups++ > MAX_N_STEP_UPS) {
712 }
713 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
714 continue;
715 }
716 /* calculate derivatives and exit function at new point */
717 (*derivs)(dydx1, y1, x1);
718 ex1 = exit_value - y1[i_exit_value];
719 /* check for change in sign of exit function */
720 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
721 if (n_to_skip == 0)
722 break;
723 else {
724 --n_to_skip;
725 is_zero = 1;
726 }
727 }
728 /* check for end of interval */
729 if (fabs(xdiff = xf - x1) < x_accuracy) {
730 /* end of the interval */
731 for (i = 0; i < n_eq; i++)
732 y_return[i] = y1[i];
733 *x0 = x1;
734 *h_rec = h_start;
735 tfree(dydx0);
736 tfree(dydx1);
737 tfree(dydx2);
738 tfree(yscale);
739 tfree(accur);
740 if (y0 != y_return)
741 tfree(y0);
742 if (y1 != y_return)
743 tfree(y1);
744 if (y2 != y_return)
745 tfree(y2);
746 return (DIFFEQ_END_OF_INTERVAL);
747 }
748 /* copy the new solution into the old variables */
749 SWAP_PTR(dydx0, dydx1);
750 SWAP_PTR(y0, y1);
751 ex0 = ex1;
752 *x0 = x1;
753 /* adjust the step size as recommended by bs_step() */
754 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
755 /* calculate new scale factors */
756 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
757 accur, n_eq);
758 } while (1);
759 *h_rec = h_start;
760
761 /* The root has been bracketed. */
762 do {
763 /* try to take a step to the position where the zero is expected */
764 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
765 x2 = *x0;
766 /* calculate new scale factors */
767 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
768 accur, n_eq);
769 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
770 yscale, n_eq, derivs, misses))
772 /* check the exit function at the new position */
773 (*derivs)(dydx2, y2, x2);
774 ex2 = exit_value - y2[i_exit_value];
775 if (fabs(ex2) < exit_accuracy) {
776 for (i = 0; i < n_eq; i++)
777 y_return[i] = y2[i];
778 *x0 = x2;
779 tfree(dydx0);
780 tfree(dydx1);
781 tfree(dydx2);
782 tfree(yscale);
783 tfree(accur);
784 if (y0 != y_return)
785 tfree(y0);
786 if (y1 != y_return)
787 tfree(y1);
788 if (y2 != y_return)
789 tfree(y2);
790 return (DIFFEQ_ZERO_FOUND);
791 }
792 /* rebracket the root */
793 if (SIGN(ex1) == SIGN(ex2)) {
794 SWAP_PTR(y1, y2);
795 SWAP_PTR(dydx1, dydx2);
796 x1 = x2;
797 ex1 = ex2;
798 } else {
799 SWAP_PTR(y0, y2);
800 SWAP_PTR(dydx0, dydx2);
801 *x0 = x2;
802 ex0 = ex2;
803 }
804 } while (1);
805}
◆ bs_odeint3()
| long bs_odeint3 | ( | double * | yif, |
| void(* | derivs )(double *dydx, double *y, double x), | ||
| long | n_eq, | ||
| double * | accuracy, | ||
| long * | accmode, | ||
| double * | tiny, | ||
| long * | misses, | ||
| double * | x0, | ||
| double | xf, | ||
| double | x_accuracy, | ||
| double | h_start, | ||
| double | h_max, | ||
| double * | h_rec, | ||
| double(* | exit_func )(double *dydx, double *y, double x), | ||
| double | exit_accuracy ) |
Integrates a system of ordinary differential equations using the Bulirsch-Stoer method with optimized internal state management.
This function is a variant of bs_odeint() that utilizes internal static variables to manage state across multiple integration steps, potentially improving performance in scenarios requiring repeated integrations.
- Parameters
-
yif Pointer to the array of initial/final values of dependent variables. Upon successful completion, it contains the final values. derivs Function pointer to compute the derivatives. It calculates dy/dx given the current state. n_eq Number of equations or dependent variables in the system. accuracy Pointer to the array specifying the desired accuracy for each dependent variable. accmode Pointer to the array specifying the accuracy control mode for each dependent variable. Modes: - 0: Fractional accuracy per step for each variable.
- 1: Fractional accuracy globally.
- 2: Absolute accuracy per step for each variable.
- 3: Absolute accuracy globally.
tiny Pointer to the array of small values representing the lower limits of significance for each dependent variable. misses Pointer to the array that counts the number of times each variable caused a step size reset due to exceeding error tolerance. x0 Pointer to the initial value of the independent variable. It is updated to the final value after integration. xf Upper limit of the independent variable to integrate up to. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size for the integration. h_max Maximum allowed step size for the integration. h_rec Pointer to the variable where the recommended step size for continuation will be stored. exit_func Function pointer to the exit condition function. It returns a value that, when zero, signals the integration to stop. exit_accuracy Desired accuracy for the exit condition function to evaluate to zero.
- Returns
- Returns a non-negative value on failure (with specific codes) or a positive value on success.
Definition at line 836 of file bsODEp.c.
857 {
858 static double *yscale;
859 static double *dydx0, *y1, *dydx1, *dydx2, *y2, *accur, *y0;
860 static long last_neq = 0;
861 double ex0, ex1, ex2, x1, x2;
862 double h_used, h_next, xdiff;
863 long i, n_step_ups = 0;
864#define MAX_N_STEP_UPS 10
865
866 if (*x0 > xf)
867 return (DIFFEQ_XI_GT_XF);
868 if (FABS(*x0 - xf) < x_accuracy)
869 return (DIFFEQ_SOLVED_ALREADY);
870
871 /* Meaning of accmode:
872 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
873 * each step for ith variable. tiny[i] is lower limit
874 * of significance for the ith variable.
875 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
876 * satisfied globally, not locally.
877 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
878 * step for the ith variable. tiny[i] is ignored.
879 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
880 * be satisfied globally, not locally.
881 */
882 for (i = 0; i < n_eq; i++) {
883 if (accmode[i] < 0 || accmode[i] > 3)
885 if (accmode[i] < 2 && tiny[i] < TINY)
886 tiny[i] = TINY;
887 misses[i] = 0;
888 }
889
890 if (last_neq < n_eq) {
891 if (last_neq != 0) {
892 tfree(y0);
893 tfree(dydx0);
894 tfree(y1);
895 tfree(dydx1);
896 tfree(y2);
897 tfree(dydx2);
898 tfree(yscale);
899 tfree(accur);
900 }
909 last_neq = n_eq;
910 }
911
912 for (i = 0; i < n_eq; i++)
913 y0[i] = yif[i];
914
915 /* calculate derivatives and exit function at the initial point */
916 (*derivs)(dydx0, y0, *x0);
917
918 /* set the scales for evaluating accuracy. yscale[i] is the
919 * absolute level of accuracy required of the next integration step
920 */
921 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
922 accuracy, accur, *x0, xf, n_eq);
923
924 ex0 = (*exit_func)(dydx0, y0, *x0);
925
926 do {
927 /* check for zero of exit function */
928 if (FABS(ex0) < exit_accuracy) {
929 for (i = 0; i < n_eq; i++)
930 yif[i] = y0[i];
931 *h_rec = h_start;
932 return (DIFFEQ_ZERO_FOUND);
933 }
934
935 /* adjust step size to stay within interval */
936 if ((xdiff = xf - *x0) < h_start)
937 h_start = xdiff;
938 /* take a step */
939 x1 = *x0;
940 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
941 yscale, n_eq, derivs, misses)) {
942 if (n_step_ups++ > MAX_N_STEP_UPS)
944 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
945 continue;
946 }
947 /* calculate derivatives and exit function at new point */
948 (*derivs)(dydx1, y1, x1);
949 ex1 = (*exit_func)(dydx1, y1, x1);
950 if (SIGN(ex0) != SIGN(ex1))
951 break;
952 /* check for end of interval */
953 if (FABS(xdiff = xf - x1) < x_accuracy) {
954 /* end of the interval */
955 for (i = 0; i < n_eq; i++)
956 yif[i] = y1[i];
957 *x0 = x1;
958 *h_rec = h_start;
959 return (DIFFEQ_END_OF_INTERVAL);
960 }
961 /* copy the new solution into the old variables */
962 SWAP_PTR(dydx0, dydx1);
963 SWAP_PTR(y0, y1);
964 ex0 = ex1;
965 *x0 = x1;
966 /* adjust the step size as recommended by bs_step() */
967 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
968 /* calculate new scale factors */
969 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
970 accur, n_eq);
971 } while (1);
972 *h_rec = h_start;
973
974 if (!exit_func) {
975 printf("failure in bs_odeint3(): solution stepped outside interval\n");
976 return (DIFFEQ_OUTSIDE_INTERVAL);
977 }
978
979 if (FABS(ex1) < exit_accuracy) {
980 for (i = 0; i < n_eq; i++)
981 yif[i] = y1[i];
982 *x0 = x1;
983 return (DIFFEQ_ZERO_FOUND);
984 }
985
986 /* The root has been bracketed. */
987 do {
988 /* try to take a step to the position where the zero is expected */
989 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
990 x2 = *x0;
991 /* calculate new scale factors */
992 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
993 accur, n_eq);
994 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
995 yscale, n_eq, derivs, misses))
997 /* check the exit function at the new position */
998 (*derivs)(dydx2, y2, x2);
999 ex2 = (*exit_func)(dydx2, y2, x2);
1000 if (FABS(ex2) < exit_accuracy) {
1001 for (i = 0; i < n_eq; i++)
1002 yif[i] = y2[i];
1003 *x0 = x2;
1004 return (DIFFEQ_ZERO_FOUND);
1005 }
1006 /* rebracket the root */
1007 if (SIGN(ex1) == SIGN(ex2)) {
1008 SWAP_PTR(y1, y2);
1009 SWAP_PTR(dydx1, dydx2);
1010 x1 = x2;
1011 ex1 = ex2;
1012 } else {
1013 SWAP_PTR(y0, y2);
1014 SWAP_PTR(dydx0, dydx2);
1015 *x0 = x2;
1016 ex0 = ex2;
1017 }
1018 } while (1);
1019}
◆ bs_odeint4()
| long bs_odeint4 | ( | double * | y0, |
| void(* | derivs )(double *dydx, double *y, double x), | ||
| long | n_eq, | ||
| double * | accuracy, | ||
| long * | accmode, | ||
| double * | tiny, | ||
| long * | misses, | ||
| double * | x0, | ||
| double | xf, | ||
| double | x_accuracy, | ||
| double | h_start, | ||
| double | h_max, | ||
| double * | h_rec, | ||
| double | exit_value, | ||
| long | i_exit_value, | ||
| double | exit_accuracy, | ||
| long | n_to_skip, | ||
| void(* | store_data )(double *dydx, double *y, double x, double exf) ) |
Integrates a system of ordinary differential equations until a specified component reaches a target value or the upper limit is met, with intermediate data storage.
This function extends bs_odeint2() by allowing the storage of intermediate integration points. It integrates the ODE system until either the independent variable reaches the specified upper limit or a particular dependent variable reaches a target value within a specified accuracy.
- Parameters
-
y0 Pointer to the array of initial values of dependent variables. Upon successful completion, it contains the final values. derivs Function pointer to compute the derivatives. It calculates dy/dx given the current state. n_eq Number of equations or dependent variables in the system. accuracy Pointer to the array specifying the desired accuracy for each dependent variable. accmode Pointer to the array specifying the accuracy control mode for each dependent variable. Modes: - 0: Fractional accuracy per step for each variable.
- 1: Fractional accuracy globally.
- 2: Absolute accuracy per step for each variable.
- 3: Absolute accuracy globally.
tiny Pointer to the array of small values representing the lower limits of significance for each dependent variable. misses Pointer to the array that counts the number of times each variable caused a step size reset due to exceeding error tolerance. x0 Pointer to the initial value of the independent variable. It is updated to the final value after integration. xf Upper limit of the independent variable to integrate up to. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size for the integration. h_max Maximum allowed step size for the integration. h_rec Pointer to the variable where the recommended step size for continuation will be stored. exit_value Target value that the specified component of the solution should reach. i_exit_value Index of the dependent variable component that is being monitored for reaching the target value. exit_accuracy Desired accuracy for the target value condition. n_to_skip Number of times the target condition can be met before integration stops. store_data Function pointer to store intermediate integration points. It is called with the current derivatives, dependent variables, independent variable, and exit function value.
- Returns
- Returns 1 on successful integration, or a non-negative value on failure.
Definition at line 1051 of file bsODEp.c.
1075 {
1076 double *y_return, *accur;
1077 double *dydx0, *y1, *dydx1, *dydx2, *y2;
1078 double ex0, ex1, ex2, x1, x2, *yscale;
1079 double h_used, h_next, xdiff;
1080 long i, n_step_ups = 0, is_zero;
1081#define MAX_N_STEP_UPS 10
1082
1083 if (*x0 > xf)
1084 return (DIFFEQ_XI_GT_XF);
1085 if (fabs(*x0 - xf) < x_accuracy)
1086 return (DIFFEQ_SOLVED_ALREADY);
1087 if (i_exit_value < 0 || i_exit_value >= n_eq)
1089
1090 /* Meaning of accmode:
1091 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
1092 * each step for ith variable. tiny[i] is lower limit
1093 * of significance for the ith variable.
1094 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
1095 * satisfied globally, not locally.
1096 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
1097 * step for the ith variable. tiny[i] is ignored.
1098 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
1099 * be satisfied globally, not locally.
1100 */
1101 for (i = 0; i < n_eq; i++) {
1102 if (accmode[i] < 0 || accmode[i] > 3)
1104 if (accmode[i] < 2 && tiny[i] < TINY)
1105 tiny[i] = TINY;
1106 misses[i] = 0;
1107 }
1108
1109 y_return = y0;
1116
1117 /* calculate derivatives and exit function at the initial point */
1118 (*derivs)(dydx0, y0, *x0);
1119
1120 /* set the scales for evaluating accuracy. yscale[i] is the
1121 * absolute level of accuracy required of the next integration step
1122 */
1124 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1125 accuracy, accur, *x0, xf, n_eq);
1126
1127 ex0 = exit_value - y0[i_exit_value];
1128 if (store_data)
1129 (*store_data)(dydx0, y0, *x0, ex0);
1130 is_zero = 0;
1131 do {
1132 /* check for zero of exit function */
1133 if (fabs(ex0) < exit_accuracy) {
1134 if (!is_zero) {
1135 if (n_to_skip == 0) {
1136 if (store_data)
1137 (*store_data)(dydx0, y0, *x0, ex0);
1138 for (i = 0; i < n_eq; i++)
1139 y_return[i] = y0[i];
1140 *h_rec = h_start;
1141 tfree(dydx0);
1142 tfree(dydx1);
1143 tfree(dydx2);
1144 tfree(yscale);
1145 tfree(accur);
1146 if (y0 != y_return)
1147 tfree(y0);
1148 if (y1 != y_return)
1149 tfree(y1);
1150 if (y2 != y_return)
1151 tfree(y2);
1152 return (DIFFEQ_ZERO_FOUND);
1153 } else {
1154 is_zero = 1;
1155 --n_to_skip;
1156 }
1157 }
1158 } else
1159 is_zero = 0;
1160 /* adjust step size to stay within interval */
1161 if ((xdiff = xf - *x0) < h_start)
1162 h_start = xdiff;
1163 /* take a step */
1164 x1 = *x0;
1165 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
1166 yscale, n_eq, derivs, misses)) {
1167 if (n_step_ups++ > MAX_N_STEP_UPS) {
1169 }
1170 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
1171 continue;
1172 }
1173 /* calculate derivatives and exit function at new point */
1174 (*derivs)(dydx1, y1, x1);
1175 ex1 = exit_value - y1[i_exit_value];
1176 if (store_data)
1177 (*store_data)(dydx1, y1, x1, ex1);
1178 /* check for change in sign of exit function */
1179 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
1180 if (n_to_skip == 0)
1181 break;
1182 else {
1183 --n_to_skip;
1184 is_zero = 1;
1185 }
1186 }
1187 /* check for end of interval */
1188 if (fabs(xdiff = xf - x1) < x_accuracy) {
1189 /* end of the interval */
1190 if (store_data) {
1191 (*derivs)(dydx1, y1, x1);
1192 ex1 = exit_value - y0[i_exit_value];
1193 (*store_data)(dydx1, y1, x1, ex1);
1194 }
1195 for (i = 0; i < n_eq; i++)
1196 y_return[i] = y1[i];
1197 *x0 = x1;
1198 *h_rec = h_start;
1199 tfree(dydx0);
1200 tfree(dydx1);
1201 tfree(dydx2);
1202 tfree(yscale);
1203 tfree(accur);
1204 if (y0 != y_return)
1205 tfree(y0);
1206 if (y1 != y_return)
1207 tfree(y1);
1208 if (y2 != y_return)
1209 tfree(y2);
1210 return (DIFFEQ_END_OF_INTERVAL);
1211 }
1212 /* copy the new solution into the old variables */
1213 SWAP_PTR(dydx0, dydx1);
1214 SWAP_PTR(y0, y1);
1215 ex0 = ex1;
1216 *x0 = x1;
1217 /* adjust the step size as recommended by bs_step() */
1218 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
1219 /* calculate new scale factors */
1220 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1221 accur, n_eq);
1222 } while (1);
1223 *h_rec = h_start;
1224
1225 /* The root has been bracketed. */
1226 do {
1227 /* try to take a step to the position where the zero is expected */
1228 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
1229 x2 = *x0;
1230 /* calculate new scale factors */
1231 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1232 accur, n_eq);
1233 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1234 yscale, n_eq, derivs, misses))
1236 /* check the exit function at the new position */
1237 (*derivs)(dydx2, y2, x2);
1238 ex2 = exit_value - y2[i_exit_value];
1239 if (fabs(ex2) < exit_accuracy) {
1240 for (i = 0; i < n_eq; i++)
1241 y_return[i] = y2[i];
1242 *x0 = x2;
1243 tfree(dydx0);
1244 tfree(dydx1);
1245 tfree(dydx2);
1246 tfree(yscale);
1247 tfree(accur);
1248 if (y0 != y_return)
1249 tfree(y0);
1250 if (y1 != y_return)
1251 tfree(y1);
1252 if (y2 != y_return)
1253 tfree(y2);
1254 return (DIFFEQ_ZERO_FOUND);
1255 }
1256 /* rebracket the root */
1257 if (SIGN(ex1) == SIGN(ex2)) {
1258 SWAP_PTR(y1, y2);
1259 SWAP_PTR(dydx1, dydx2);
1260 x1 = x2;
1261 ex1 = ex2;
1262 } else {
1263 SWAP_PTR(y0, y2);
1264 SWAP_PTR(dydx0, dydx2);
1265 *x0 = x2;
1266 ex0 = ex2;
1267 }
1268 } while (1);
1269}
◆ bs_qctune()
| void bs_qctune | ( | double | newStepIncreaseFactor, |
| double | newStepDecreaseFactor ) |
◆ bs_step()
| long bs_step | ( | double * | yFinal, |
| double * | x, | ||
| double * | yInitial, | ||
| double * | dydxInitial, | ||
| double | step, | ||
| double * | stepUsed, | ||
| double * | stepRecommended, | ||
| double * | yScale, | ||
| long | equations, | ||
| void(* | derivs )(double *dydx, double *y, double x), | ||
| long * | misses ) |
Definition at line 47 of file bsODEp.c.
60 {
61 static double **solution, *hSqr, *yLast, *yError;
62 long i, j, iWorst = 0, code, nuse;
63 double maxError, error, yInterp;
64 static long mmidSteps[IMAX] = {2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96};
65 static long lastEquations = 0;
66
67 if (equations > lastEquations) {
68 if (lastEquations != 0) {
69 free(yLast);
70 free(yError);
71 free(hSqr);
73 }
78 lastEquations = equations;
79 }
80
81 do {
82#if DEBUG
83 printf("step = %e\n", step);
84#endif
85 for (i = 0; i < IMAX; i++) {
86 mmid(yInitial, dydxInitial, equations, *x, step, mmidSteps[i], yFinal, derivs);
87 hSqr[i % NUSE] = sqr(step / mmidSteps[i]);
88 nuse = i > NUSE ? NUSE : i;
89 for (j = 0; j < equations; j++) {
90 /* store the jth component of the new solution, possibly throwing out the oldest solution */
91 solution[j][i % NUSE] = yFinal[j];
92 if (nuse > 1)
93 /* Interpolate the solution value at h^2 = 0 */
94 yInterp = LagrangeInterp(hSqr, solution[j], nuse, 0.0, &code);
95 else
96 /* just copy modified midpoint value */
97 yInterp = yFinal[j];
98 if (i)
99 /* compute difference between new solution and that from last step */
100 yError[j] = yInterp - yLast[j];
101 /* save the new solution for the next iteration */
102 yLast[j] = yInterp;
103#if DEBUG
104 printf("i = %ld, j = %ld: yFinal = %.10e, yError = %.10e, yScale = %.10e\n",
105 i, j, yFinal[j], yInterp, yError[j], yScale[j]);
106#endif
107 }
108 if (i) {
109 maxError = 0.0;
110 for (j = 0; j < equations; j++)
111 if (maxError < (error = FABS(yError[j] / yScale[j]))) {
112 iWorst = j;
113 maxError = error;
114 }
115#if DEBUG
116 printf("maxError = %e, iWorst = %ld\n", maxError, iWorst);
117#endif
118 if (maxError < 1.0) {
119 *x += step;
120 *stepRecommended = *stepUsed = step;
121 if (i == NUSE - 1)
122 /* had a hard time, so recommend smaller step */
123 *stepRecommended *= stepDecreaseFactor;
124 else {
125 /* increase the step to make better use of extrapolation */
126 *stepRecommended *= stepIncreaseFactor / sqrt(maxError);
127 }
128#if DEBUG
129 printf("returning with i=%ld, stepUsed=%e, stepRec=%e\n", i, *stepUsed, *stepRecommended);
130#endif
131 for (j = 0; j < equations; j++)
132 yFinal[j] = yLast[j];
133 return (1);
134 }
135 misses[iWorst]++;
136 }
137 }
138
139 /* method failed, so reduce the step of the modified-midpoint integrations */
140 step *= 0.25;
141 for (i = 0; i < (IMAX - NUSE) / 2; i++)
142 step /= 2.0;
143 } while ((*x + step) != *x);
144 fprintf(stderr, "error: step size underflow in bs_step()--step reduced to %e\n", step);
145 return 0;
146}
void ** array_2d(uint64_t size, uint64_t lower1, uint64_t upper1, uint64_t lower2, uint64_t upper2)
Allocates a 2D array with specified lower and upper indices for both dimensions.
Definition array.c:275
int free_array_2d(void **array, uint64_t size, uint64_t lower1, uint64_t upper1, uint64_t lower2, uint64_t upper2)
Frees a 2D array and its associated memory.
Definition array.c:327
double LagrangeInterp(double *x, double *f, long order1, double x0, long *returnCode)
Performs Lagrange interpolation of data.
Definition interp.c:115
◆ initial_scale_factors_dp()
| void initial_scale_factors_dp | ( | double * | yscale, |
| double * | y0, | ||
| double * | dydx0, | ||
| double | h_start, | ||
| double * | tiny, | ||
| long * | accmode, | ||
| double * | accuracy, | ||
| double * | accur, | ||
| double | x0, | ||
| double | xf, | ||
| long | n_eq ) |
Definition at line 70 of file rkODE.c.
80 {
81 int i;
82
83 for (i = 0; i < n_eq; i++) {
84 if ((accur[i] = accuracy[i]) <= 0) {
85 printf("error: accuracy[%d] = %e (initial_scale_factors_dp)\n", i, accuracy[i]);
86 abort();
87 }
88 switch (accmode[i]) {
89 case -1: /* no accuracy control */
90 yscale[i] = DBL_MAX;
91 break;
92 case 0: /* fractional local accuracy specified */
93 yscale[i] = (y0[i] + dydx0[i] * h_start + tiny[i]) * accur[i];
94 break;
95 case 1: /* fractional global accuracy specified */
96 yscale[i] = (dydx0[i] * h_start + tiny[i]) * accur[i];
97 break;
98 case 2: /* absolute local accuracy specified */
99 yscale[i] = accur[i];
100 break;
101 case 3: /* absolute global accuracy specified */
102 yscale[i] = (accur[i] /= (xf - x0)) * h_start;
103 break;
104 default:
105 printf("error: accmode[%d] = %ld (initial_scale_factors_dp)\n", i, accmode[i]);
106 abort();
107 break;
108 }
109 if (yscale[i] <= 0) {
110 printf("error: yscale[%d] = %e (initial_scale_factors_dp)\n", i, yscale[i]);
111 abort();
112 }
113 }
114}
◆ new_scale_factors_dp()
| void new_scale_factors_dp | ( | double * | yscale, |
| double * | y0, | ||
| double * | dydx0, | ||
| double | h_start, | ||
| double * | tiny, | ||
| long * | accmode, | ||
| double * | accuracy, | ||
| long | n_eq ) |
Definition at line 28 of file rkODE.c.
36 {
37 int i;
38
39 for (i = 0; i < n_eq; i++) {
40 switch (accmode[i]) {
41 case -1: /* no accuracy control */
42 yscale[i] = DBL_MAX;
43 break;
44 case 0: /* fractional local accuracy specified */
45 yscale[i] =
46 (y0[i] + dydx0[i] * h_start + tiny[i]) * accuracy[i];
47 break;
48 case 1: /* fractional global accuracy specified */
49 yscale[i] =
50 (dydx0[i] * h_start + tiny[i]) * accuracy[i];
51 break;
52 case 2: /* absolute local accuracy specified */
53 yscale[i] = accuracy[i];
54 break;
55 case 3: /* absolute global accuracy specified */
56 yscale[i] = accuracy[i] * h_start;
57 break;
58 default:
59 printf("error: accmode[%d] = %ld (new_scale_factors_dp)\n", i, accmode[i]);
60 abort();
61 break;
62 }
63 if (yscale[i] <= 0) {
64 printf("error: yscale[%d] = %e\n", i, yscale[i]);
65 abort();
66 }
67 }
68}
Variable Documentation
◆ stepDecreaseFactor
◆ stepIncreaseFactor
Generated by
1.12.0