bsODEp.c File Reference
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.
#include "mdb.h"#include "mdb_thread.h"Go to the source code of this file.
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)(double *yp, 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) |
| 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. | |
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 189 of file bsODEp.c.
213 {
214 double *y_return, *accur;
215 double *dydx0, *y1, *dydx1, *dydx2, *y2;
216 double ex0, ex1, ex2, x1, x2, *yscale;
217 double h_used, h_next, xdiff;
218 long i, n_step_ups = 0, is_zero;
219#define MAX_N_STEP_UPS 10
220
221 if (*x0 > xf)
222 return (DIFFEQ_XI_GT_XF);
223 if (FABS(*x0 - xf) < x_accuracy)
224 return (DIFFEQ_SOLVED_ALREADY);
225
226 /* Meaning of accmode:
227 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
228 * each step for ith variable. tiny[i] is lower limit
229 * of significance for the ith variable.
230 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
231 * satisfied globally, not locally.
232 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
233 * step for the ith variable. tiny[i] is ignored.
234 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
235 * be satisfied globally, not locally.
236 */
237 for (i = 0; i < n_eq; i++) {
238 if (accmode[i] < 0 || accmode[i] > 3)
240 if (accmode[i] < 2 && tiny[i] < TINY)
241 tiny[i] = TINY;
242 misses[i] = 0;
243 }
244
245 y_return = y0;
252
253 /* calculate derivatives and exit function at the initial point */
254 (*derivs)(dydx0, y0, *x0);
255
256 /* set the scales for evaluating accuracy. yscale[i] is the
257 * absolute level of accuracy required of the next integration step
258 */
260 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
261 accuracy, accur, *x0, xf, n_eq);
262
263 ex0 = exit_func ? (*exit_func)(dydx0, y0, *x0) : 0;
264 if (store_data)
265 (*store_data)(dydx0, y0, *x0, ex0);
266 is_zero = 0;
267
268 do {
269 /* check for zero of exit function */
270 if (exit_func && FABS(ex0) < exit_accuracy) {
271 if (!is_zero) {
272 if (n_to_skip == 0) {
273 if (store_data)
274 (*store_data)(dydx0, y0, *x0, ex0);
275 for (i = 0; i < n_eq; i++)
276 y_return[i] = y0[i];
277 *h_rec = h_start;
278 tfree(dydx0);
279 tfree(dydx1);
280 tfree(dydx2);
281 tfree(yscale);
282 tfree(accur);
283 if (y0 != y_return)
284 tfree(y0);
285 if (y1 != y_return)
286 tfree(y1);
287 if (y2 != y_return)
288 tfree(y2);
289 return (DIFFEQ_ZERO_FOUND);
290 } else {
291 is_zero = 1;
292 --n_to_skip;
293 }
294 }
295 } else
296 is_zero = 0;
297 /* adjust step size to stay within interval */
298 if ((xdiff = xf - *x0) < h_start)
299 h_start = xdiff;
300 /* take a step */
301 x1 = *x0;
302 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
303 yscale, n_eq, derivs, misses)) {
304 if (n_step_ups++ > MAX_N_STEP_UPS)
306 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
307 continue;
308 }
309 /* calculate derivatives and exit function at new point */
310 (*derivs)(dydx1, y1, x1);
311 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
312 if (store_data)
313 (*store_data)(dydx1, y1, x1, ex1);
314 /* check for change in sign of exit function */
315 if (exit_func && SIGN(ex0) != SIGN(ex1) && !is_zero) {
316 if (n_to_skip == 0)
317 break;
318 else {
319 --n_to_skip;
320 is_zero = 1;
321 }
322 }
323 /* check for end of interval */
324 if (FABS(xdiff = xf - x1) < x_accuracy) {
325 /* end of the interval */
326 if (store_data) {
327 (*derivs)(dydx1, y1, x1);
328 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
329 (*store_data)(dydx1, y1, x1, ex1);
330 }
331 for (i = 0; i < n_eq; i++)
332 y_return[i] = y1[i];
333 *x0 = x1;
334 *h_rec = h_start;
335 tfree(dydx0);
336 tfree(dydx1);
337 tfree(dydx2);
338 tfree(yscale);
339 tfree(accur);
340 if (y0 != y_return)
341 tfree(y0);
342 if (y1 != y_return)
343 tfree(y1);
344 if (y2 != y_return)
345 tfree(y2);
346 return (DIFFEQ_END_OF_INTERVAL);
347 }
348 /* copy the new solution into the old variables */
349 SWAP_PTR(dydx0, dydx1);
350 SWAP_PTR(y0, y1);
351 ex0 = ex1;
352 *x0 = x1;
353 /* adjust the step size as recommended by bs_step() */
354 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
355 /* calculate new scale factors */
356 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
357 accur, n_eq);
358 } while (1);
359 *h_rec = h_start;
360
361 if (!exit_func) {
362 printf("failure in bs_odeint(): solution stepped outside interval\n");
363 tfree(dydx0);
364 tfree(dydx1);
365 tfree(dydx2);
366 tfree(yscale);
367 tfree(accur);
368 if (y0 != y_return)
369 tfree(y0);
370 if (y1 != y_return)
371 tfree(y1);
372 if (y2 != y_return)
373 tfree(y2);
374 return (DIFFEQ_OUTSIDE_INTERVAL);
375 }
376
377 if (FABS(ex1) < exit_accuracy) {
378 for (i = 0; i < n_eq; i++)
379 y_return[i] = y1[i];
380 *x0 = x1;
381 tfree(dydx0);
382 tfree(dydx1);
383 tfree(dydx2);
384 tfree(yscale);
385 tfree(accur);
386 if (y0 != y_return)
387 tfree(y0);
388 if (y1 != y_return)
389 tfree(y1);
390 if (y2 != y_return)
391 tfree(y2);
392 return (DIFFEQ_ZERO_FOUND);
393 }
394
395 /* The root has been bracketed. */
396 do {
397 /* try to take a step to the position where the zero is expected */
398 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
399 x2 = *x0;
400 /* calculate new scale factors */
401 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
402 accur, n_eq);
403 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
404 yscale, n_eq, derivs, misses))
406 /* check the exit function at the new position */
407 (*derivs)(dydx2, y2, x2);
408 ex2 = (*exit_func)(dydx2, y2, x2);
409 if (FABS(ex2) < exit_accuracy) {
410 for (i = 0; i < n_eq; i++)
411 y_return[i] = y2[i];
412 *x0 = x2;
413 tfree(dydx0);
414 tfree(dydx1);
415 tfree(dydx2);
416 tfree(yscale);
417 tfree(accur);
418 if (y0 != y_return)
419 tfree(y0);
420 if (y1 != y_return)
421 tfree(y1);
422 if (y2 != y_return)
423 tfree(y2);
424 return (DIFFEQ_ZERO_FOUND);
425 }
426 /* rebracket the root */
427 if (SIGN(ex1) == SIGN(ex2)) {
428 SWAP_PTR(y1, y2);
429 SWAP_PTR(dydx1, dydx2);
430 x1 = x2;
431 ex1 = ex2;
432 } else {
433 SWAP_PTR(y0, y2);
434 SWAP_PTR(dydx0, dydx2);
435 *x0 = x2;
436 ex0 = ex2;
437 }
438 } while (1);
439}
int tfree(void *ptr)
Frees a memory block and records the deallocation if tracking is enabled.
Definition array.c:243
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:65
◆ bs_odeint1()
| long bs_odeint1 | ( | double * | y0, |
| void(* | derivs )(double *yp, 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 ) |
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 466 of file bsODEp.c.
483 {
484 double *y_return;
485 double *dydx0, *y1, *dydx1, *dydx2, *y2;
486 double x1, *yscale, *accur;
487 double h_used, h_next, xdiff;
488 long i, n_step_ups = 0;
489#define MAX_N_STEP_UPS 10
490
491 if (*x0 > xf)
492 return (DIFFEQ_XI_GT_XF);
493 if (fabs(*x0 - xf) < x_accuracy)
494 return (DIFFEQ_SOLVED_ALREADY);
495
496 /* Meaning of accmode:
497 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
498 * each step for ith variable. tiny[i] is lower limit
499 * of significance for the ith variable.
500 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
501 * satisfied globally, not locally.
502 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
503 * step for the ith variable. tiny[i] is ignored.
504 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
505 * be satisfied globally, not locally.
506 */
507 for (i = 0; i < n_eq; i++) {
508 if (accmode[i] < 0 || accmode[i] > 3)
510 if (accmode[i] < 2 && tiny[i] < TINY)
511 tiny[i] = TINY;
512 misses[i] = 0;
513 }
514
515 y_return = y0;
522
523 /* calculate derivatives at the initial point */
524 (*derivs)(dydx0, y0, *x0);
525
526 /* set the scales for evaluating accuracy. yscale[i] is the
527 * absolute level of accuracy required of the next integration step
528 */
530 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
531 accuracy, accur, *x0, xf, n_eq);
532
533 do {
534 /* adjust step size to stay within interval */
535 if ((xdiff = xf - *x0) < h_start)
536 h_start = xdiff;
537 /* take a step */
538 x1 = *x0;
539 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
540 yscale, n_eq, derivs, misses)) {
541 if (n_step_ups++ > MAX_N_STEP_UPS)
543 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
544 continue;
545 }
546 /* check for end of interval */
547 if (fabs(xdiff = xf - x1) < x_accuracy) {
548 /* end of the interval */
549 for (i = 0; i < n_eq; i++)
550 y_return[i] = y1[i];
551 *x0 = x1;
552 *h_rec = h_start;
553 tfree(dydx0);
554 tfree(dydx1);
555 tfree(dydx2);
556 tfree(yscale);
557 tfree(accur);
558 if (y0 != y_return)
559 tfree(y0);
560 if (y1 != y_return)
561 tfree(y1);
562 if (y2 != y_return)
563 tfree(y2);
564 return (DIFFEQ_END_OF_INTERVAL);
565 }
566 /* calculate derivatives at new point */
567 (*derivs)(dydx1, y1, x1);
568 /* copy the new solution into the old variables */
569 SWAP_PTR(dydx0, dydx1);
570 SWAP_PTR(y0, y1);
571 *x0 = x1;
572 /* adjust the step size as recommended by bs_step() */
573 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
574 /* calculate new scale factors */
575 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
576 accur, n_eq);
577 } while (1);
578}
◆ 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 609 of file bsODEp.c.
632 {
633 double *y_return, *accur;
634 double *dydx0, *y1, *dydx1, *dydx2, *y2;
635 double ex0, ex1, ex2, x1, x2, *yscale;
636 double h_used, h_next, xdiff;
637 long i, n_step_ups = 0, is_zero;
638#define MAX_N_STEP_UPS 10
639
640 if (*x0 > xf)
641 return (DIFFEQ_XI_GT_XF);
642 if (fabs(*x0 - xf) < x_accuracy)
643 return (DIFFEQ_SOLVED_ALREADY);
644 if (i_exit_value < 0 || i_exit_value >= n_eq)
646
647 /* Meaning of accmode:
648 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
649 * each step for ith variable. tiny[i] is lower limit
650 * of significance for the ith variable.
651 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
652 * satisfied globally, not locally.
653 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
654 * step for the ith variable. tiny[i] is ignored.
655 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
656 * be satisfied globally, not locally.
657 */
658 for (i = 0; i < n_eq; i++) {
659 if (accmode[i] < 0 || accmode[i] > 3)
661 if (accmode[i] < 2 && tiny[i] < TINY)
662 tiny[i] = TINY;
663 misses[i] = 0;
664 }
665
666 y_return = y0;
673
674 /* calculate derivatives and exit function at the initial point */
675 (*derivs)(dydx0, y0, *x0);
676
677 /* set the scales for evaluating accuracy. yscale[i] is the
678 * absolute level of accuracy required of the next integration step
679 */
681 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
682 accuracy, accur, *x0, xf, n_eq);
683
684 ex0 = exit_value - y0[i_exit_value];
685 is_zero = 0;
686 do {
687 /* check for zero of exit function */
688 if (fabs(ex0) < exit_accuracy) {
689 if (!is_zero) {
690 if (n_to_skip == 0) {
691 for (i = 0; i < n_eq; i++)
692 y_return[i] = y0[i];
693 *h_rec = h_start;
694 tfree(dydx0);
695 tfree(dydx1);
696 tfree(dydx2);
697 tfree(yscale);
698 tfree(accur);
699 if (y0 != y_return)
700 tfree(y0);
701 if (y1 != y_return)
702 tfree(y1);
703 if (y2 != y_return)
704 tfree(y2);
705 return (DIFFEQ_ZERO_FOUND);
706 } else {
707 is_zero = 1;
708 --n_to_skip;
709 }
710 }
711 } else
712 is_zero = 0;
713 /* adjust step size to stay within interval */
714 if ((xdiff = xf - *x0) < h_start)
715 h_start = xdiff;
716 /* take a step */
717 x1 = *x0;
718 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
719 yscale, n_eq, derivs, misses)) {
720 if (n_step_ups++ > MAX_N_STEP_UPS) {
722 }
723 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
724 continue;
725 }
726 /* calculate derivatives and exit function at new point */
727 (*derivs)(dydx1, y1, x1);
728 ex1 = exit_value - y1[i_exit_value];
729 /* check for change in sign of exit function */
730 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
731 if (n_to_skip == 0)
732 break;
733 else {
734 --n_to_skip;
735 is_zero = 1;
736 }
737 }
738 /* check for end of interval */
739 if (fabs(xdiff = xf - x1) < x_accuracy) {
740 /* end of the interval */
741 for (i = 0; i < n_eq; i++)
742 y_return[i] = y1[i];
743 *x0 = x1;
744 *h_rec = h_start;
745 tfree(dydx0);
746 tfree(dydx1);
747 tfree(dydx2);
748 tfree(yscale);
749 tfree(accur);
750 if (y0 != y_return)
751 tfree(y0);
752 if (y1 != y_return)
753 tfree(y1);
754 if (y2 != y_return)
755 tfree(y2);
756 return (DIFFEQ_END_OF_INTERVAL);
757 }
758 /* copy the new solution into the old variables */
759 SWAP_PTR(dydx0, dydx1);
760 SWAP_PTR(y0, y1);
761 ex0 = ex1;
762 *x0 = x1;
763 /* adjust the step size as recommended by bs_step() */
764 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
765 /* calculate new scale factors */
766 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
767 accur, n_eq);
768 } while (1);
769 *h_rec = h_start;
770
771 /* The root has been bracketed. */
772 do {
773 /* try to take a step to the position where the zero is expected */
774 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
775 x2 = *x0;
776 /* calculate new scale factors */
777 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
778 accur, n_eq);
779 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
780 yscale, n_eq, derivs, misses))
782 /* check the exit function at the new position */
783 (*derivs)(dydx2, y2, x2);
784 ex2 = exit_value - y2[i_exit_value];
785 if (fabs(ex2) < exit_accuracy) {
786 for (i = 0; i < n_eq; i++)
787 y_return[i] = y2[i];
788 *x0 = x2;
789 tfree(dydx0);
790 tfree(dydx1);
791 tfree(dydx2);
792 tfree(yscale);
793 tfree(accur);
794 if (y0 != y_return)
795 tfree(y0);
796 if (y1 != y_return)
797 tfree(y1);
798 if (y2 != y_return)
799 tfree(y2);
800 return (DIFFEQ_ZERO_FOUND);
801 }
802 /* rebracket the root */
803 if (SIGN(ex1) == SIGN(ex2)) {
804 SWAP_PTR(y1, y2);
805 SWAP_PTR(dydx1, dydx2);
806 x1 = x2;
807 ex1 = ex2;
808 } else {
809 SWAP_PTR(y0, y2);
810 SWAP_PTR(dydx0, dydx2);
811 *x0 = x2;
812 ex0 = ex2;
813 }
814 } while (1);
815}
◆ 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 846 of file bsODEp.c.
867 {
868 static MDB_THREAD_LOCAL double *yscale = NULL;
869 static MDB_THREAD_LOCAL double *dydx0 = NULL, *y1 = NULL, *dydx1 = NULL, *dydx2 = NULL, *y2 = NULL, *accur = NULL, *y0 = NULL;
870 static MDB_THREAD_LOCAL long last_neq = 0;
871 double ex0, ex1, ex2, x1, x2;
872 double h_used, h_next, xdiff;
873 long i, n_step_ups = 0;
874#define MAX_N_STEP_UPS 10
875
876 if (*x0 > xf)
877 return (DIFFEQ_XI_GT_XF);
878 if (FABS(*x0 - xf) < x_accuracy)
879 return (DIFFEQ_SOLVED_ALREADY);
880
881 /* Meaning of accmode:
882 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
883 * each step for ith variable. tiny[i] is lower limit
884 * of significance for the ith variable.
885 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
886 * satisfied globally, not locally.
887 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
888 * step for the ith variable. tiny[i] is ignored.
889 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
890 * be satisfied globally, not locally.
891 */
892 for (i = 0; i < n_eq; i++) {
893 if (accmode[i] < 0 || accmode[i] > 3)
895 if (accmode[i] < 2 && tiny[i] < TINY)
896 tiny[i] = TINY;
897 misses[i] = 0;
898 }
899
900 if (last_neq < n_eq) {
901 if (last_neq != 0) {
902 tfree(y0);
903 tfree(dydx0);
904 tfree(y1);
905 tfree(dydx1);
906 tfree(y2);
907 tfree(dydx2);
908 tfree(yscale);
909 tfree(accur);
910 }
919 last_neq = n_eq;
920 }
921
922 for (i = 0; i < n_eq; i++)
923 y0[i] = yif[i];
924
925 /* calculate derivatives and exit function at the initial point */
926 (*derivs)(dydx0, y0, *x0);
927
928 /* set the scales for evaluating accuracy. yscale[i] is the
929 * absolute level of accuracy required of the next integration step
930 */
931 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
932 accuracy, accur, *x0, xf, n_eq);
933
934 ex0 = (*exit_func)(dydx0, y0, *x0);
935
936 do {
937 /* check for zero of exit function */
938 if (FABS(ex0) < exit_accuracy) {
939 for (i = 0; i < n_eq; i++)
940 yif[i] = y0[i];
941 *h_rec = h_start;
942 return (DIFFEQ_ZERO_FOUND);
943 }
944
945 /* adjust step size to stay within interval */
946 if ((xdiff = xf - *x0) < h_start)
947 h_start = xdiff;
948 /* take a step */
949 x1 = *x0;
950 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
951 yscale, n_eq, derivs, misses)) {
952 if (n_step_ups++ > MAX_N_STEP_UPS)
954 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
955 continue;
956 }
957 /* calculate derivatives and exit function at new point */
958 (*derivs)(dydx1, y1, x1);
959 ex1 = (*exit_func)(dydx1, y1, x1);
960 if (SIGN(ex0) != SIGN(ex1))
961 break;
962 /* check for end of interval */
963 if (FABS(xdiff = xf - x1) < x_accuracy) {
964 /* end of the interval */
965 for (i = 0; i < n_eq; i++)
966 yif[i] = y1[i];
967 *x0 = x1;
968 *h_rec = h_start;
969 return (DIFFEQ_END_OF_INTERVAL);
970 }
971 /* copy the new solution into the old variables */
972 SWAP_PTR(dydx0, dydx1);
973 SWAP_PTR(y0, y1);
974 ex0 = ex1;
975 *x0 = x1;
976 /* adjust the step size as recommended by bs_step() */
977 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
978 /* calculate new scale factors */
979 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
980 accur, n_eq);
981 } while (1);
982 *h_rec = h_start;
983
984 if (!exit_func) {
985 printf("failure in bs_odeint3(): solution stepped outside interval\n");
986 return (DIFFEQ_OUTSIDE_INTERVAL);
987 }
988
989 if (FABS(ex1) < exit_accuracy) {
990 for (i = 0; i < n_eq; i++)
991 yif[i] = y1[i];
992 *x0 = x1;
993 return (DIFFEQ_ZERO_FOUND);
994 }
995
996 /* The root has been bracketed. */
997 do {
998 /* try to take a step to the position where the zero is expected */
999 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
1000 x2 = *x0;
1001 /* calculate new scale factors */
1002 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1003 accur, n_eq);
1004 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1005 yscale, n_eq, derivs, misses))
1007 /* check the exit function at the new position */
1008 (*derivs)(dydx2, y2, x2);
1009 ex2 = (*exit_func)(dydx2, y2, x2);
1010 if (FABS(ex2) < exit_accuracy) {
1011 for (i = 0; i < n_eq; i++)
1012 yif[i] = y2[i];
1013 *x0 = x2;
1014 return (DIFFEQ_ZERO_FOUND);
1015 }
1016 /* rebracket the root */
1017 if (SIGN(ex1) == SIGN(ex2)) {
1018 SWAP_PTR(y1, y2);
1019 SWAP_PTR(dydx1, dydx2);
1020 x1 = x2;
1021 ex1 = ex2;
1022 } else {
1023 SWAP_PTR(y0, y2);
1024 SWAP_PTR(dydx0, dydx2);
1025 *x0 = x2;
1026 ex0 = ex2;
1027 }
1028 } while (1);
1029}
◆ 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 1061 of file bsODEp.c.
1085 {
1086 double *y_return, *accur;
1087 double *dydx0, *y1, *dydx1, *dydx2, *y2;
1088 double ex0, ex1, ex2, x1, x2, *yscale;
1089 double h_used, h_next, xdiff;
1090 long i, n_step_ups = 0, is_zero;
1091#define MAX_N_STEP_UPS 10
1092
1093 if (*x0 > xf)
1094 return (DIFFEQ_XI_GT_XF);
1095 if (fabs(*x0 - xf) < x_accuracy)
1096 return (DIFFEQ_SOLVED_ALREADY);
1097 if (i_exit_value < 0 || i_exit_value >= n_eq)
1099
1100 /* Meaning of accmode:
1101 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
1102 * each step for ith variable. tiny[i] is lower limit
1103 * of significance for the ith variable.
1104 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
1105 * satisfied globally, not locally.
1106 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
1107 * step for the ith variable. tiny[i] is ignored.
1108 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
1109 * be satisfied globally, not locally.
1110 */
1111 for (i = 0; i < n_eq; i++) {
1112 if (accmode[i] < 0 || accmode[i] > 3)
1114 if (accmode[i] < 2 && tiny[i] < TINY)
1115 tiny[i] = TINY;
1116 misses[i] = 0;
1117 }
1118
1119 y_return = y0;
1126
1127 /* calculate derivatives and exit function at the initial point */
1128 (*derivs)(dydx0, y0, *x0);
1129
1130 /* set the scales for evaluating accuracy. yscale[i] is the
1131 * absolute level of accuracy required of the next integration step
1132 */
1134 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1135 accuracy, accur, *x0, xf, n_eq);
1136
1137 ex0 = exit_value - y0[i_exit_value];
1138 if (store_data)
1139 (*store_data)(dydx0, y0, *x0, ex0);
1140 is_zero = 0;
1141 do {
1142 /* check for zero of exit function */
1143 if (fabs(ex0) < exit_accuracy) {
1144 if (!is_zero) {
1145 if (n_to_skip == 0) {
1146 if (store_data)
1147 (*store_data)(dydx0, y0, *x0, ex0);
1148 for (i = 0; i < n_eq; i++)
1149 y_return[i] = y0[i];
1150 *h_rec = h_start;
1151 tfree(dydx0);
1152 tfree(dydx1);
1153 tfree(dydx2);
1154 tfree(yscale);
1155 tfree(accur);
1156 if (y0 != y_return)
1157 tfree(y0);
1158 if (y1 != y_return)
1159 tfree(y1);
1160 if (y2 != y_return)
1161 tfree(y2);
1162 return (DIFFEQ_ZERO_FOUND);
1163 } else {
1164 is_zero = 1;
1165 --n_to_skip;
1166 }
1167 }
1168 } else
1169 is_zero = 0;
1170 /* adjust step size to stay within interval */
1171 if ((xdiff = xf - *x0) < h_start)
1172 h_start = xdiff;
1173 /* take a step */
1174 x1 = *x0;
1175 if (!bs_step(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
1176 yscale, n_eq, derivs, misses)) {
1177 if (n_step_ups++ > MAX_N_STEP_UPS) {
1179 }
1180 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
1181 continue;
1182 }
1183 /* calculate derivatives and exit function at new point */
1184 (*derivs)(dydx1, y1, x1);
1185 ex1 = exit_value - y1[i_exit_value];
1186 if (store_data)
1187 (*store_data)(dydx1, y1, x1, ex1);
1188 /* check for change in sign of exit function */
1189 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
1190 if (n_to_skip == 0)
1191 break;
1192 else {
1193 --n_to_skip;
1194 is_zero = 1;
1195 }
1196 }
1197 /* check for end of interval */
1198 if (fabs(xdiff = xf - x1) < x_accuracy) {
1199 /* end of the interval */
1200 if (store_data) {
1201 (*derivs)(dydx1, y1, x1);
1202 ex1 = exit_value - y0[i_exit_value];
1203 (*store_data)(dydx1, y1, x1, ex1);
1204 }
1205 for (i = 0; i < n_eq; i++)
1206 y_return[i] = y1[i];
1207 *x0 = x1;
1208 *h_rec = h_start;
1209 tfree(dydx0);
1210 tfree(dydx1);
1211 tfree(dydx2);
1212 tfree(yscale);
1213 tfree(accur);
1214 if (y0 != y_return)
1215 tfree(y0);
1216 if (y1 != y_return)
1217 tfree(y1);
1218 if (y2 != y_return)
1219 tfree(y2);
1220 return (DIFFEQ_END_OF_INTERVAL);
1221 }
1222 /* copy the new solution into the old variables */
1223 SWAP_PTR(dydx0, dydx1);
1224 SWAP_PTR(y0, y1);
1225 ex0 = ex1;
1226 *x0 = x1;
1227 /* adjust the step size as recommended by bs_step() */
1228 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
1229 /* calculate new scale factors */
1230 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1231 accur, n_eq);
1232 } while (1);
1233 *h_rec = h_start;
1234
1235 /* The root has been bracketed. */
1236 do {
1237 /* try to take a step to the position where the zero is expected */
1238 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0 + TINY);
1239 x2 = *x0;
1240 /* calculate new scale factors */
1241 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1242 accur, n_eq);
1243 if (!bs_step(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1244 yscale, n_eq, derivs, misses))
1246 /* check the exit function at the new position */
1247 (*derivs)(dydx2, y2, x2);
1248 ex2 = exit_value - y2[i_exit_value];
1249 if (fabs(ex2) < exit_accuracy) {
1250 for (i = 0; i < n_eq; i++)
1251 y_return[i] = y2[i];
1252 *x0 = x2;
1253 tfree(dydx0);
1254 tfree(dydx1);
1255 tfree(dydx2);
1256 tfree(yscale);
1257 tfree(accur);
1258 if (y0 != y_return)
1259 tfree(y0);
1260 if (y1 != y_return)
1261 tfree(y1);
1262 if (y2 != y_return)
1263 tfree(y2);
1264 return (DIFFEQ_ZERO_FOUND);
1265 }
1266 /* rebracket the root */
1267 if (SIGN(ex1) == SIGN(ex2)) {
1268 SWAP_PTR(y1, y2);
1269 SWAP_PTR(dydx1, dydx2);
1270 x1 = x2;
1271 ex1 = ex2;
1272 } else {
1273 SWAP_PTR(y0, y2);
1274 SWAP_PTR(dydx0, dydx2);
1275 *x0 = x2;
1276 ex0 = ex2;
1277 }
1278 } while (1);
1279}
◆ 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 51 of file bsODEp.c.
64 {
65 static MDB_THREAD_LOCAL double **solution = NULL, *hSqr = NULL, *yLast = NULL, *yError = NULL;
66 long i, j, iWorst = 0, code, nuse;
67 double maxError, error, yInterp;
68 double localStepIncreaseFactor, localStepDecreaseFactor;
69 static const long mmidSteps[IMAX] = {2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96};
70 static MDB_THREAD_LOCAL long lastEquations = 0;
71
72 mdb_thread_lock(&bs_qctune_lock);
73 localStepIncreaseFactor = stepIncreaseFactor;
74 localStepDecreaseFactor = stepDecreaseFactor;
75 mdb_thread_unlock(&bs_qctune_lock);
76
77 if (equations > lastEquations) {
78 if (lastEquations != 0) {
79 free(yLast);
80 free(yError);
81 free(hSqr);
83 }
88 lastEquations = equations;
89 }
90
91 do {
92#if DEBUG
93 printf("step = %e\n", step);
94#endif
95 for (i = 0; i < IMAX; i++) {
96 mmid(yInitial, dydxInitial, equations, *x, step, mmidSteps[i], yFinal, derivs);
97 hSqr[i % NUSE] = sqr(step / mmidSteps[i]);
98 nuse = i > NUSE ? NUSE : i;
99 for (j = 0; j < equations; j++) {
100 /* store the jth component of the new solution, possibly throwing out the oldest solution */
101 solution[j][i % NUSE] = yFinal[j];
102 if (nuse > 1)
103 /* Interpolate the solution value at h^2 = 0 */
104 yInterp = LagrangeInterp(hSqr, solution[j], nuse, 0.0, &code);
105 else
106 /* just copy modified midpoint value */
107 yInterp = yFinal[j];
108 if (i)
109 /* compute difference between new solution and that from last step */
110 yError[j] = yInterp - yLast[j];
111 /* save the new solution for the next iteration */
112 yLast[j] = yInterp;
113#if DEBUG
114 printf("i = %ld, j = %ld: yFinal = %.10e, yError = %.10e, yScale = %.10e\n",
115 i, j, yFinal[j], yInterp, yError[j], yScale[j]);
116#endif
117 }
118 if (i) {
119 maxError = 0.0;
120 for (j = 0; j < equations; j++)
121 if (maxError < (error = FABS(yError[j] / yScale[j]))) {
122 iWorst = j;
123 maxError = error;
124 }
125#if DEBUG
126 printf("maxError = %e, iWorst = %ld\n", maxError, iWorst);
127#endif
128 if (maxError < 1.0) {
129 *x += step;
130 *stepRecommended = *stepUsed = step;
131 if (i == NUSE - 1)
132 /* had a hard time, so recommend smaller step */
133 *stepRecommended *= localStepDecreaseFactor;
134 else {
135 /* increase the step to make better use of extrapolation */
136 *stepRecommended *= localStepIncreaseFactor / sqrt(maxError);
137 }
138#if DEBUG
139 printf("returning with i=%ld, stepUsed=%e, stepRec=%e\n", i, *stepUsed, *stepRecommended);
140#endif
141 for (j = 0; j < equations; j++)
142 yFinal[j] = yLast[j];
143 return (1);
144 }
145 misses[iWorst]++;
146 }
147 }
148
149 /* method failed, so reduce the step of the modified-midpoint integrations */
150 step *= 0.25;
151 for (i = 0; i < (IMAX - NUSE) / 2; i++)
152 step /= 2.0;
153 } while ((*x + step) != *x);
154 fprintf(stderr, "error: step size underflow in bs_step()--step reduced to %e\n", step);
155 return 0;
156}
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:290
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:342
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 71 of file rkODE.c.
81 {
82 int i;
83
84 for (i = 0; i < n_eq; i++) {
85 if ((accur[i] = accuracy[i]) <= 0) {
86 printf("error: accuracy[%d] = %e (initial_scale_factors_dp)\n", i, accuracy[i]);
87 abort();
88 }
89 switch (accmode[i]) {
90 case -1: /* no accuracy control */
91 yscale[i] = DBL_MAX;
92 break;
93 case 0: /* fractional local accuracy specified */
94 yscale[i] = (y0[i] + dydx0[i] * h_start + tiny[i]) * accur[i];
95 break;
96 case 1: /* fractional global accuracy specified */
97 yscale[i] = (dydx0[i] * h_start + tiny[i]) * accur[i];
98 break;
99 case 2: /* absolute local accuracy specified */
100 yscale[i] = accur[i];
101 break;
102 case 3: /* absolute global accuracy specified */
103 yscale[i] = (accur[i] /= (xf - x0)) * h_start;
104 break;
105 default:
106 printf("error: accmode[%d] = %ld (initial_scale_factors_dp)\n", i, accmode[i]);
107 abort();
108 break;
109 }
110 if (yscale[i] <= 0) {
111 printf("error: yscale[%d] = %e (initial_scale_factors_dp)\n", i, yscale[i]);
112 abort();
113 }
114 }
115}
◆ 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 29 of file rkODE.c.
37 {
38 int i;
39
40 for (i = 0; i < n_eq; i++) {
41 switch (accmode[i]) {
42 case -1: /* no accuracy control */
43 yscale[i] = DBL_MAX;
44 break;
45 case 0: /* fractional local accuracy specified */
46 yscale[i] =
47 (y0[i] + dydx0[i] * h_start + tiny[i]) * accuracy[i];
48 break;
49 case 1: /* fractional global accuracy specified */
50 yscale[i] =
51 (dydx0[i] * h_start + tiny[i]) * accuracy[i];
52 break;
53 case 2: /* absolute local accuracy specified */
54 yscale[i] = accuracy[i];
55 break;
56 case 3: /* absolute global accuracy specified */
57 yscale[i] = accuracy[i] * h_start;
58 break;
59 default:
60 printf("error: accmode[%d] = %ld (new_scale_factors_dp)\n", i, accmode[i]);
61 abort();
62 break;
63 }
64 if (yscale[i] <= 0) {
65 printf("error: yscale[%d] = %e\n", i, yscale[i]);
66 abort();
67 }
68 }
69}
Generated by
1.12.0