rkODE.c File Reference
Detailed Description
Fourth-order Runge-Kutta ODE integration routines (double-precision version).
This file provides functions for integrating ordinary differential equations using fourth-order Runge-Kutta methods. It includes adaptive step-size control and supports various integration scenarios.
- 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 rkODE.c.
#include "mdb.h"#include "mdb_thread.h"#include <float.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 | report_state_dp (FILE *fp, double *y, double *dydx, double *yscale, long *misses, double x, double h, long n_eq) |
| void | rk4_step (double *yf, double x, double *yi, double *dydx, double h, long n_eq, void(*derivs)(double *dydx, double *y, double x)) |
| void | rk4_qctune (double newSafetyMargin, double newIncreasePower, double newDecreasePower, double newMaxIncreaseFactor) |
| long | rk_qcstep (double *yFinal, double *x, double *yInitial, double *dydxInitial, double hInput, double *hUsed, double *hRecommended, double *yScale, long equations, void(*derivs)(double *dydx, double *y, double x), long *misses) |
| long | rk_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 exval)) |
| Integrate a set of ODEs until the upper limit or an exit condition is met. | |
| long | rk_odeint1 (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) |
| Integrate ODEs without exit conditions or intermediate output. | |
| long | rk_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) |
| Integrate ODEs until a specific component reaches a target value. | |
| long | rk_odeint_na (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, double h_max, double *h_rec) |
| Integrate ODEs without adaptive step-size control. | |
| long | rk_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) |
| Integrate ODEs with an exit condition. | |
| long | rk_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)) |
| Integrate ODEs until a specific component reaches a target value with intermediate storage. | |
| long | rk_odeint3_na (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_step, double h_max, double *h_rec, double(*exit_func)(double *dydx, double *y, double x), double exit_accuracy, void(*stochastic)(double *y, double x, double h)) |
| Integrate ODEs without adaptive step-size and with stochastic processes. | |
Function Documentation
◆ 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}
◆ report_state_dp()
| void report_state_dp | ( | FILE * | fp, |
| double * | y, | ||
| double * | dydx, | ||
| double * | yscale, | ||
| long * | misses, | ||
| double | x, | ||
| double | h, | ||
| long | n_eq ) |
Definition at line 117 of file rkODE.c.
118 {
119 int i;
120 fputs("integration state:\n", fp);
121 fprintf(fp, "%ld equations, indep.var.=%e, step size=%e",
122 n_eq, x, h);
123 fprintf(fp, "\ny : ");
124 for (i = 0; i < n_eq; i++)
125 fprintf(fp, "%e ", y[i]);
126 fprintf(fp, "\ndydx : ");
127 for (i = 0; i < n_eq; i++)
128 fprintf(fp, "%e ", dydx[i]);
129 fprintf(fp, "\ntol.scale: ");
130 for (i = 0; i < n_eq; i++)
131 fprintf(fp, "%e ", yscale[i]);
132 fprintf(fp, "\nmisses : ");
133 for (i = 0; i < n_eq; i++)
134 fprintf(fp, "%ld ", misses[i]);
135}
◆ rk4_qctune()
| void rk4_qctune | ( | double | newSafetyMargin, |
| double | newIncreasePower, | ||
| double | newDecreasePower, | ||
| double | newMaxIncreaseFactor ) |
Definition at line 222 of file rkODE.c.
223 {
224 mdb_thread_lock(&rk_qctune_lock);
225 if (newSafetyMargin > 0 && newSafetyMargin < 1)
226 safetyMargin = newSafetyMargin;
227 if (newIncreasePower > 0)
228 increasePower = newIncreasePower;
229 if (newDecreasePower > 0)
230 decreasePower = newDecreasePower;
231 if (newMaxIncreaseFactor > 1)
232 maxIncreaseFactor = newMaxIncreaseFactor;
233 mdb_thread_unlock(&rk_qctune_lock);
234}
◆ rk4_step()
| void rk4_step | ( | double * | yf, |
| double | x, | ||
| double * | yi, | ||
| double * | dydx, | ||
| double | h, | ||
| long | n_eq, | ||
| void(* | derivs )(double *dydx, double *y, double x) ) |
Definition at line 142 of file rkODE.c.
151 {
152 static MDB_THREAD_LOCAL long last_n_eq = 0;
153 static MDB_THREAD_LOCAL double *k1 = NULL, *k2 = NULL, *k3 = NULL, *yTemp = NULL, *dydxTemp = NULL;
154 double x1;
155 long i;
156
157 /* for speed, I avoid reallocating the arrays unless it is necessary */
158 if (last_n_eq < n_eq) {
159 if (last_n_eq != 0) {
160 free(k1);
161 free(k2);
162 free(k3);
163 free(yTemp);
164 free(dydxTemp);
165 }
166 last_n_eq = n_eq;
172 }
173
174 /* yf = yi + k1/6 + k2/3 + k3/3 + k4/6 where
175 k1 = h*dydt(x , y)
176 k2 = h*dydt(x+h/2, y+k1/2)
177 k3 = h*dydt(x+h/2, y+k2/2)
178 k4 = h*dydt(x+h , y+k3)
179 */
180
181 /* compute k1 and yTemp for k2 computation */
182 for (i = 0; i < n_eq; i++) {
183 k1[i] = h * dydx[i];
184 yTemp[i] = yi[i] + k1[i] / 2;
185 }
186
187 /* compute k2 and yTemp for k3 computation */
188 x1 = x + h / 2;
189 (*derivs)(dydxTemp, yTemp, x1);
190 for (i = 0; i < n_eq; i++) {
191 k2[i] = h * dydxTemp[i];
192 yTemp[i] = yi[i] + k2[i] / 2;
193 }
194
195 /* compute k3 and yTemp for k4 computation */
196 (*derivs)(dydxTemp, yTemp, x1);
197 for (i = 0; i < n_eq; i++) {
198 k3[i] = h * dydxTemp[i];
199 yTemp[i] = yi[i] + k3[i];
200 }
201
202 /* compute k4 and put results into yf array */
203 x1 = x + h;
204 (*derivs)(dydxTemp, yTemp, x1);
205 for (i = 0; i < n_eq; i++)
206 yf[i] = yi[i] + (k1[i] / 2 + k2[i] + k3[i] + h * dydxTemp[i] / 2) / 3;
207}
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:65
◆ rk_odeint()
| long rk_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 exval) ) |
Integrate a set of ODEs until the upper limit or an exit condition is met.
Integrates a system of ordinary differential equations using adaptive step-size control until the independent variable reaches the upper limit or a user-defined exit condition is satisfied.
- Parameters
-
y0 Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Desired accuracies for each variable. accmode Accuracy control modes for each variable. tiny Small values relative to what's important for each variable. misses Array tracking the number of step size reductions per variable. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size. h_max Maximum step size allowed. h_rec Pointer to store the recommended step size for continuation. exit_func Function to determine when to stop integration based on a condition. exit_accuracy Accuracy required for the exit condition. n_to_skip Number of zeros of the exit function to skip before returning. store_data Function to store intermediate data points.
- Returns
- DIFFEQ_ZERO_FOUND (>=1) on success, or error code (<=0) on failure.
Definition at line 409 of file rkODE.c.
433 {
434 double *y_return;
435 double *dydx0, *y1, *dydx1, *dydx2, *y2;
436 double ex0, ex1, ex2, *yscale, *accur;
437 double h_used, h_next, x1, x2, xdiff;
438 long i, n_exit_iterations, n_step_ups = 0, is_zero;
439#define MAX_N_STEP_UPS 10
440
441 if (*x0 > xf)
442 return (DIFFEQ_XI_GT_XF);
443 if (FABS(*x0 - xf) < x_accuracy)
444 return (DIFFEQ_SOLVED_ALREADY);
445
446 /* Meaning of accmode:
447 * accmode = -1 -> no accuracy control
448 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
449 * each step for ith variable. tiny[i] is lower limit
450 * of significance for the ith variable.
451 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
452 * satisfied globally, not locally.
453 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
454 * step for the ith variable. tiny[i] is ignored.
455 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
456 * be satisfied globally, not locally.
457 */
458 for (i = 0; i < n_eq; i++) {
459 if (accmode[i] < -1 || accmode[i] > 3)
461 if (accmode[i] < 2 && tiny[i] < TINY)
462 tiny[i] = TINY;
463 misses[i] = 0;
464 }
465
466 y_return = y0;
473
474 /* calculate derivatives and exit function at the initial point */
475 (*derivs)(dydx0, y0, *x0);
476
477 /* set the scales for evaluating accuracy. yscale[i] is the
478 * absolute level of accuracy required of the next integration step.
479 */
481 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
482 accuracy, accur, *x0, xf, n_eq);
483
484 ex0 = exit_func ? (*exit_func)(dydx0, y0, *x0) : 0;
485 if (store_data)
486 (*store_data)(dydx0, y0, *x0, ex0);
487 is_zero = 0;
488
489 do {
490 /* check for zero of exit function */
491 if (exit_func && FABS(ex0) < exit_accuracy) {
492 if (!is_zero) {
493 if (n_to_skip == 0) {
494 if (store_data)
495 (*store_data)(dydx0, y0, *x0, ex0);
496 for (i = 0; i < n_eq; i++)
497 y_return[i] = y0[i];
498 *h_rec = h_start;
499 tfree(dydx0);
500 tfree(dydx1);
501 tfree(dydx2);
502 tfree(yscale);
503 tfree(accur);
504 if (y0 != y_return)
505 tfree(y0);
506 if (y1 != y_return)
507 tfree(y1);
508 if (y2 != y_return)
509 tfree(y2);
510 return (DIFFEQ_ZERO_FOUND);
511 } else {
512 is_zero = 1;
513 --n_to_skip;
514 }
515 }
516 } else
517 is_zero = 0;
518 /* adjust step size to stay within interval */
519 if ((xdiff = xf - *x0) < h_start)
520 h_start = xdiff;
521 /* take a step */
522 x1 = *x0;
523 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
524 yscale, n_eq, derivs, misses)) {
525 if (n_step_ups++ > MAX_N_STEP_UPS)
527 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
528 continue;
529 }
530 /* calculate derivatives and exit function at new point */
531 (*derivs)(dydx1, y1, x1);
532 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
533 if (store_data)
534 (*store_data)(dydx1, y1, x1, ex1);
535 /* check for change in sign of exit function */
536 if (exit_func && SIGN(ex0) != SIGN(ex1) && !is_zero) {
537 if (n_to_skip == 0)
538 break;
539 else {
540 --n_to_skip;
541 is_zero = 1;
542 }
543 }
544
545 /* check for end of interval */
546 if (FABS(xdiff = xf - x1) < x_accuracy) {
547 /* end of the interval */
548 if (store_data) {
549 (*derivs)(dydx1, y1, x1);
550 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
551 (*store_data)(dydx1, y1, x1, ex1);
552 }
553 for (i = 0; i < n_eq; i++)
554 y_return[i] = y1[i];
555 *x0 = x1;
556 *h_rec = h_start;
557 tfree(dydx0);
558 tfree(dydx1);
559 tfree(dydx2);
560 tfree(yscale);
561 tfree(accur);
562 if (y0 != y_return)
563 tfree(y0);
564 if (y1 != y_return)
565 tfree(y1);
566 if (y2 != y_return)
567 tfree(y2);
568 return (DIFFEQ_END_OF_INTERVAL);
569 }
570
571 /* copy the new solution into the old variables */
572 SWAP_PTR(dydx0, dydx1);
573 SWAP_PTR(y0, y1);
574 ex0 = ex1;
575 *x0 = x1;
576 /* adjust the step size as recommended by rk_qcstep() */
577 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
578 /* calculate new scale factors */
579 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
580 accur, n_eq);
581
582 } while (1);
583 *h_rec = h_start;
584
585 if (!exit_func) {
586 printf("failure in rk_odeint(): solution stepped outside interval\n");
587 tfree(dydx0);
588 tfree(dydx1);
589 tfree(dydx2);
590 tfree(yscale);
591 tfree(accur);
592 if (y0 != y_return)
593 tfree(y0);
594 if (y1 != y_return)
595 tfree(y1);
596 if (y2 != y_return)
597 tfree(y2);
598 return (DIFFEQ_OUTSIDE_INTERVAL);
599 }
600
601 /* The root has been bracketed. */
602 n_exit_iterations = MAX_EXIT_ITERATIONS;
603 do {
604 /* try to take a step to the position where the zero is expected */
605 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
606 x2 = *x0;
607 /* calculate new scale factors */
608 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
609 accur, n_eq);
610 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
611 yscale, n_eq, derivs, misses))
613 /* check the exit function at the new position */
614 (*derivs)(dydx2, y2, x2);
615 ex2 = (*exit_func)(dydx2, y2, x2);
616 if (FABS(ex2) < exit_accuracy) {
617 for (i = 0; i < n_eq; i++)
618 y_return[i] = y2[i];
619 *x0 = x2;
620 tfree(dydx0);
621 tfree(dydx1);
622 tfree(dydx2);
623 tfree(yscale);
624 tfree(accur);
625 if (y0 != y_return)
626 tfree(y0);
627 if (y1 != y_return)
628 tfree(y1);
629 if (y2 != y_return)
630 tfree(y2);
631 return (DIFFEQ_ZERO_FOUND);
632 }
633 /* rebracket the root */
634 if (SIGN(ex1) == SIGN(ex2)) {
635 SWAP_PTR(y1, y2);
636 SWAP_PTR(dydx1, dydx2);
637 x1 = x2;
638 ex1 = ex2;
639 } else {
640 SWAP_PTR(y0, y2);
641 SWAP_PTR(dydx0, dydx2);
642 *x0 = x2;
643 ex0 = ex2;
644 }
645 } while (n_exit_iterations--);
646 return (DIFFEQ_EXIT_COND_FAILED);
647}
int tfree(void *ptr)
Frees a memory block and records the deallocation if tracking is enabled.
Definition array.c:243
◆ rk_odeint1()
| long rk_odeint1 | ( | 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 ) |
Integrate ODEs without exit conditions or intermediate output.
Integrates a system of ordinary differential equations until the upper limit of the independent variable is reached. This function does not monitor exit conditions or store intermediate data, making it faster for simple integrations.
- Parameters
-
y0 Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Ignored in this function. accmode Ignored in this function. tiny Ignored in this function. misses Ignored in this function. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Ignored in this function. h_start Step size. h_max Maximum step size allowed. h_rec Pointer to store the recommended step size for continuation.
- Returns
- DIFFEQ_END_OF_INTERVAL (1) on success, or an error code (<=0) on failure.
Definition at line 671 of file rkODE.c.
689 {
690 double *y_return;
691 double *dydx0, *y1, *dydx1, *dydx2, *y2;
692 double *yscale, *accur;
693 double x1;
694 double h_used, h_next, xdiff;
695 long i, n_step_ups = 0;
696#define MAX_N_STEP_UPS 10
697
698 if (*x0 > xf)
699 return (DIFFEQ_XI_GT_XF);
700 if (FABS(*x0 - xf) < x_accuracy)
701 return (DIFFEQ_SOLVED_ALREADY);
702
703 /* Meaning of accmode:
704 * accmode = -1 -> no accuracy control
705 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
706 * each step for ith variable. tiny[i] is lower limit
707 * of significance for the ith variable.
708 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
709 * satisfied globally, not locally.
710 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
711 * step for the ith variable. tiny[i] is ignored.
712 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
713 * be satisfied globally, not locally.
714 */
715 for (i = 0; i < n_eq; i++) {
716 if (accmode[i] < -1 || accmode[i] > 3)
718 if (accmode[i] < 2 && tiny[i] < TINY)
719 tiny[i] = TINY;
720 misses[i] = 0;
721 }
722
723 y_return = y0;
730
731 /* calculate derivatives at the initial point */
732 (*derivs)(dydx0, y0, *x0);
733
734 /* set the scales for evaluating accuracy. yscale[i] is the
735 * absolute level of accuracy required of the next integration step
736 */
738 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
739 accuracy, accur, *x0, xf, n_eq);
740
741 do {
742 /* adjust step size to stay within interval */
743 if ((xdiff = xf - *x0) < h_start)
744 h_start = xdiff;
745 /* take a step */
746 x1 = *x0;
747 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
748 yscale, n_eq, derivs, misses)) {
749 if (n_step_ups++ > MAX_N_STEP_UPS) {
750 puts("error: cannot take step (rk_odeint1--1)");
751 printf("xf = %.16e x0 = %.16e\n", xf, *x0);
752 printf("h_start = %.16e h_used = %.16e\n", h_start, h_used);
753 puts("dump of integration state:");
754 puts(" variable value derivative scale misses");
755 puts("---------------------------------------------------------------");
756 for (i = 0; i < n_eq; i++)
757 printf(" %5ld %13.6e %13.6e %13.6e %5ld \n",
758 i, y0[i], dydx0[i], yscale[i], misses[i]);
759 exit(1);
760 }
761 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
762 continue;
763 }
764 /* check for end of interval */
765 if (FABS(xdiff = xf - x1) < x_accuracy) {
766 /* end of the interval */
767 for (i = 0; i < n_eq; i++)
768 y_return[i] = y1[i];
769 *x0 = x1;
770 *h_rec = h_start;
771 tfree(dydx0);
772 tfree(dydx1);
773 tfree(dydx2);
774 tfree(yscale);
775 tfree(accur);
776 if (y0 != y_return)
777 tfree(y0);
778 if (y1 != y_return)
779 tfree(y1);
780 if (y2 != y_return)
781 tfree(y2);
782 return (DIFFEQ_END_OF_INTERVAL);
783 }
784 /* calculate derivatives at new point */
785 (*derivs)(dydx1, y1, x1);
786 /* copy the new solution into the old variables */
787 SWAP_PTR(dydx0, dydx1);
788 SWAP_PTR(y0, y1);
789 *x0 = x1;
790 /* adjust the step size as recommended by rk_qcstep() */
791 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
792 /* calculate new scale factors */
793 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
794 accur, n_eq);
795 } while (1);
796}
◆ rk_odeint2()
| long rk_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 ) |
Integrate ODEs until a specific component reaches a target value.
Integrates a system of ordinary differential equations until a specified component of the solution reaches a target value within a given accuracy. This function does not monitor general exit conditions or store intermediate data.
- Parameters
-
y0 Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Desired accuracies for each variable. accmode Accuracy control modes for each variable. tiny Small values relative to what's important for each variable. misses Array tracking the number of step size reductions per variable. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size. h_max Maximum step size allowed. h_rec Pointer to store the recommended step size for continuation. exit_value Target value for the specified component. i_exit_value Index of the component to monitor. exit_accuracy Accuracy required for the target value. n_to_skip Number of target value crossings to skip before stopping.
- Returns
- DIFFEQ_ZERO_FOUND (1) on success, or an error code (<=0) on failure.
Definition at line 824 of file rkODE.c.
847 {
848 double *y_return, *accur;
849 double *dydx0, *y1, *dydx1, *dydx2, *y2;
850 double ex0, ex1, ex2, x1, x2, *yscale;
851 double h_used, h_next, xdiff;
852 long i, n_exit_iterations, n_step_ups = 0, is_zero;
853#define MAX_N_STEP_UPS 10
854
855 if (*x0 > xf)
856 return (DIFFEQ_XI_GT_XF);
857 if (FABS(*x0 - xf) < x_accuracy)
858 return (DIFFEQ_SOLVED_ALREADY);
859 if (i_exit_value < 0 || i_exit_value >= n_eq)
861
862 /* Meaning of accmode:
863 * accmode = -1 -> no accuracy control
864 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
865 * each step for ith variable. tiny[i] is lower limit
866 * of significance for the ith variable.
867 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
868 * satisfied globally, not locally.
869 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
870 * step for the ith variable. tiny[i] is ignored.
871 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
872 * be satisfied globally, not locally.
873 */
874 for (i = 0; i < n_eq; i++) {
875 if (accmode[i] < -1 || accmode[i] > 3)
877 if (accmode[i] < 2 && tiny[i] < TINY)
878 tiny[i] = TINY;
879 misses[i] = 0;
880 }
881
882 y_return = y0;
889
890 /* calculate derivatives and exit function at the initial point */
891 (*derivs)(dydx0, y0, *x0);
892
893 /* set the scales for evaluating accuracy. yscale[i] is the
894 * absolute level of accuracy required of the next integration step
895 */
897 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
898 accuracy, accur, *x0, xf, n_eq);
899
900 ex0 = exit_value - y0[i_exit_value];
901 is_zero = 0;
902 do {
903 /* check for zero of exit function */
904 if (FABS(ex0) < exit_accuracy) {
905 if (!is_zero) {
906 if (n_to_skip == 0) {
907 for (i = 0; i < n_eq; i++)
908 y_return[i] = y0[i];
909 *h_rec = h_start;
910 tfree(dydx0);
911 tfree(dydx1);
912 tfree(dydx2);
913 tfree(yscale);
914 tfree(accur);
915 if (y0 != y_return)
916 tfree(y0);
917 if (y1 != y_return)
918 tfree(y1);
919 if (y2 != y_return)
920 tfree(y2);
921 return (DIFFEQ_ZERO_FOUND);
922 } else {
923 is_zero = 1;
924 --n_to_skip;
925 }
926 }
927 } else
928 is_zero = 0;
929 /* adjust step size to stay within interval */
930 if ((xdiff = xf - *x0) < h_start)
931 h_start = xdiff;
932 /* take a step */
933 x1 = *x0;
934 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
935 yscale, n_eq, derivs, misses)) {
936 if (n_step_ups++ > MAX_N_STEP_UPS) {
938 }
939 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
940 continue;
941 }
942 /* calculate derivatives and exit function at new point */
943 (*derivs)(dydx1, y1, x1);
944 ex1 = exit_value - y1[i_exit_value];
945 /* check for change in sign of exit function */
946 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
947 if (n_to_skip == 0)
948 break;
949 else {
950 --n_to_skip;
951 is_zero = 1;
952 }
953 }
954 /* check for end of interval */
955 if (FABS(xdiff = xf - x1) < x_accuracy) {
956 /* end of the interval */
957 for (i = 0; i < n_eq; i++)
958 y_return[i] = y1[i];
959 *x0 = x1;
960 *h_rec = h_start;
961 tfree(dydx0);
962 tfree(dydx1);
963 tfree(dydx2);
964 tfree(yscale);
965 tfree(accur);
966 if (y0 != y_return)
967 tfree(y0);
968 if (y1 != y_return)
969 tfree(y1);
970 if (y2 != y_return)
971 tfree(y2);
972 return (DIFFEQ_END_OF_INTERVAL);
973 }
974 /* copy the new solution into the old variables */
975 SWAP_PTR(dydx0, dydx1);
976 SWAP_PTR(y0, y1);
977 ex0 = ex1;
978 *x0 = x1;
979 /* adjust the step size as recommended by rk_qcstep() */
980 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
981 /* calculate new scale factors */
982 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
983 accur, n_eq);
984
985 } while (1);
986 *h_rec = h_start;
987
988 /* The root has been bracketed. */
989 n_exit_iterations = MAX_EXIT_ITERATIONS;
990 do {
991 /* try to take a step to the position where the zero is expected */
992 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
993 x2 = *x0;
994 /* calculate new scale factors */
995 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
996 accur, n_eq);
997 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
998 yscale, n_eq, derivs, misses))
1000 /* check the exit function at the new position */
1001 (*derivs)(dydx2, y2, x2);
1002 ex2 = exit_value - y2[i_exit_value];
1003 if (FABS(ex2) < exit_accuracy) {
1004 for (i = 0; i < n_eq; i++)
1005 y_return[i] = y2[i];
1006 *x0 = x2;
1007 tfree(dydx0);
1008 tfree(dydx1);
1009 tfree(dydx2);
1010 tfree(yscale);
1011 tfree(accur);
1012 if (y0 != y_return)
1013 tfree(y0);
1014 if (y1 != y_return)
1015 tfree(y1);
1016 if (y2 != y_return)
1017 tfree(y2);
1018 return (DIFFEQ_ZERO_FOUND);
1019 }
1020 /* rebracket the root */
1021 if (SIGN(ex1) == SIGN(ex2)) {
1022 SWAP_PTR(y1, y2);
1023 SWAP_PTR(dydx1, dydx2);
1024 x1 = x2;
1025 ex1 = ex2;
1026 } else {
1027 SWAP_PTR(y0, y2);
1028 SWAP_PTR(dydx0, dydx2);
1029 *x0 = x2;
1030 ex0 = ex2;
1031 }
1032 } while (n_exit_iterations--);
1033 return (DIFFEQ_EXIT_COND_FAILED);
1034}
◆ rk_odeint3()
| long rk_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 ) |
Integrate ODEs with an exit condition.
Integrates a system of ordinary differential equations until the upper limit of the independent variable is reached or a user-defined exit condition is satisfied.
- Parameters
-
yif Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Desired accuracies for each variable. accmode Accuracy control modes for each variable. tiny Small values relative to what's important for each variable. misses Array tracking the number of step size reductions per variable. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size. h_max Maximum step size allowed. h_rec Pointer to store the recommended step size for continuation. exit_func Function to determine when to stop integration based on a condition. exit_accuracy Accuracy required for the exit condition.
- Returns
- DIFFEQ_ZERO_FOUND (>=1) on success, or an error code (<=0) on failure.
Definition at line 1159 of file rkODE.c.
1180 {
1181 static MDB_THREAD_LOCAL double *y0 = NULL, *yscale = NULL;
1182 static MDB_THREAD_LOCAL double *dydx0 = NULL, *y1 = NULL, *dydx1 = NULL, *dydx2 = NULL, *y2 = NULL, *accur = NULL;
1183 static MDB_THREAD_LOCAL long last_neq = 0;
1184 double ex0, ex1, ex2, x1, x2;
1185 double h_used, h_next, xdiff;
1186 long i, n_exit_iterations, n_step_ups = 0;
1187#define MAX_N_STEP_UPS 10
1188
1189 if (*x0 > xf)
1190 return (DIFFEQ_XI_GT_XF);
1191 if (FABS(*x0 - xf) < x_accuracy)
1192 return (DIFFEQ_SOLVED_ALREADY);
1193
1194 /* Meaning of accmode:
1195 * accmode = -1 -> no accuracy control
1196 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
1197 * each step for ith variable. tiny[i] is lower limit
1198 * of significance for the ith variable.
1199 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
1200 * satisfied globally, not locally.
1201 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
1202 * step for the ith variable. tiny[i] is ignored.
1203 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
1204 * be satisfied globally, not locally.
1205 */
1206 for (i = 0; i < n_eq; i++) {
1207 if (accmode[i] < -1 || accmode[i] > 3)
1209 if (accmode[i] < 2 && tiny[i] < TINY)
1210 tiny[i] = TINY;
1211 misses[i] = 0;
1212 }
1213
1214 if (last_neq < n_eq) {
1215 if (last_neq != 0) {
1216 tfree(y0);
1217 tfree(dydx0);
1218 tfree(y1);
1219 tfree(dydx1);
1220 tfree(y2);
1221 tfree(dydx2);
1222 tfree(yscale);
1223 tfree(accur);
1224 }
1233 last_neq = n_eq;
1234 }
1235
1236 for (i = 0; i < n_eq; i++)
1237 y0[i] = yif[i];
1238
1239 /* calculate derivatives and exit function at the initial point */
1240 (*derivs)(dydx0, y0, *x0);
1241
1242 /* set the scales for evaluating accuracy. yscale[i] is the
1243 * absolute level of accuracy required of the next integration step
1244 */
1245 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1246 accuracy, accur, *x0, xf, n_eq);
1247
1248 ex0 = (*exit_func)(dydx0, y0, *x0);
1249
1250 do {
1251 /* check for zero of exit function */
1252 if (FABS(ex0) < exit_accuracy) {
1253 for (i = 0; i < n_eq; i++)
1254 yif[i] = y0[i];
1255 *h_rec = h_start;
1256 return (DIFFEQ_ZERO_FOUND);
1257 }
1258
1259 /* adjust step size to stay within interval */
1260 if ((xdiff = xf - *x0) < h_start)
1261 h_start = xdiff;
1262 /* take a step */
1263 x1 = *x0;
1264 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
1265 yscale, n_eq, derivs, misses)) {
1266 if (n_step_ups++ > MAX_N_STEP_UPS)
1268 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
1269 continue;
1270 }
1271 /* calculate derivatives and exit function at new point */
1272 (*derivs)(dydx1, y1, x1);
1273 ex1 = (*exit_func)(dydx1, y1, x1);
1274 if (SIGN(ex0) != SIGN(ex1))
1275 break;
1276 /* check for end of interval */
1277 if (FABS(xdiff = xf - x1) < x_accuracy) {
1278 /* end of the interval */
1279 for (i = 0; i < n_eq; i++)
1280 yif[i] = y1[i];
1281 *x0 = x1;
1282 *h_rec = h_start;
1283 return (DIFFEQ_END_OF_INTERVAL);
1284 }
1285 /* copy the new solution into the old variables */
1286 SWAP_PTR(dydx0, dydx1);
1287 SWAP_PTR(y0, y1);
1288 ex0 = ex1;
1289 *x0 = x1;
1290 /* adjust the step size as recommended by rk_qcstep() */
1291 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
1292 /* calculate new scale factors */
1293 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1294 accur, n_eq);
1295 } while (1);
1296 *h_rec = h_start;
1297
1298 if (!exit_func) {
1299 printf("failure in rk_odeint3(): solution stepped outside interval\n");
1300 return (DIFFEQ_OUTSIDE_INTERVAL);
1301 }
1302
1303 if (FABS(ex1) < exit_accuracy) {
1304 for (i = 0; i < n_eq; i++)
1305 yif[i] = y1[i];
1306 *x0 = x1;
1307 return (DIFFEQ_ZERO_FOUND);
1308 }
1309
1310 /* The root has been bracketed. */
1311 n_exit_iterations = MAX_EXIT_ITERATIONS;
1312 do {
1313 /* try to take a step to the position where the zero is expected */
1314 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
1315 x2 = *x0;
1316 /* calculate new scale factors */
1317 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1318 accur, n_eq);
1319 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1320 yscale, n_eq, derivs, misses))
1322 /* check the exit function at the new position */
1323 (*derivs)(dydx2, y2, x2);
1324 ex2 = (*exit_func)(dydx2, y2, x2);
1325 if (FABS(ex2) < exit_accuracy) {
1326 for (i = 0; i < n_eq; i++)
1327 yif[i] = y2[i];
1328 *x0 = x2;
1329 return (DIFFEQ_ZERO_FOUND);
1330 }
1331 /* rebracket the root */
1332 if (SIGN(ex1) == SIGN(ex2)) {
1333 SWAP_PTR(y1, y2);
1334 SWAP_PTR(dydx1, dydx2);
1335 x1 = x2;
1336 ex1 = ex2;
1337 } else {
1338 SWAP_PTR(y0, y2);
1339 SWAP_PTR(dydx0, dydx2);
1340 *x0 = x2;
1341 ex0 = ex2;
1342 }
1343 } while (n_exit_iterations--);
1344 return (DIFFEQ_EXIT_COND_FAILED);
1345}
◆ rk_odeint3_na()
| long rk_odeint3_na | ( | 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_step, | ||
| double | h_max, | ||
| double * | h_rec, | ||
| double(* | exit_func )(double *dydx, double *y, double x), | ||
| double | exit_accuracy, | ||
| void(* | stochastic )(double *y, double x, double h) ) |
Integrate ODEs without adaptive step-size and with stochastic processes.
Integrates a system of ordinary differential equations without using adaptive step-size control and allows for the inclusion of stochastic processes.
- Parameters
-
yif Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Ignored in this function. accmode Ignored in this function. tiny Ignored in this function. misses Ignored in this function. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Desired accuracy for the final value of the independent variable. h_step Step size. h_max Ignored in this function. h_rec Ignored in this function. exit_func Function to determine when to stop integration based on a condition. exit_accuracy Accuracy required for the exit condition. stochastic Function to add stochastic processes to the solution.
- Returns
- DIFFEQ_ZERO_FOUND (>=1) on success, or an error code (<=0) on failure.
Definition at line 1621 of file rkODE.c.
1642 {
1643 static MDB_THREAD_LOCAL double *y0 = NULL, *yscale = NULL;
1644 static MDB_THREAD_LOCAL double *dydx0 = NULL, *y1 = NULL, *dydx1 = NULL, *dydx2 = NULL, *y2 = NULL, *accur = NULL;
1645 static MDB_THREAD_LOCAL long last_neq = 0;
1646 double ex0, ex1, ex2, x1, x2;
1647 double xdiff;
1648 long i, n_exit_iterations;
1649#define MAX_N_STEP_UPS 10
1650
1651 if (*x0 > xf)
1652 return (DIFFEQ_XI_GT_XF);
1653 if (FABS(*x0 - xf) < x_accuracy)
1654 return (DIFFEQ_SOLVED_ALREADY);
1655
1656 if (last_neq < n_eq) {
1657 if (last_neq != 0) {
1658 tfree(y0);
1659 tfree(dydx0);
1660 tfree(y1);
1661 tfree(dydx1);
1662 tfree(y2);
1663 tfree(dydx2);
1664 tfree(yscale);
1665 tfree(accur);
1666 }
1673 last_neq = n_eq;
1674 }
1675
1676 for (i = 0; i < n_eq; i++)
1677 y0[i] = yif[i];
1678
1679 /* calculate derivatives and exit function at the initial point */
1680 (*derivs)(dydx0, y0, *x0);
1681
1682 ex0 = (*exit_func)(dydx0, y0, *x0);
1683
1684 do {
1685 /* check for zero of exit function */
1686 if (FABS(ex0) < exit_accuracy) {
1687 for (i = 0; i < n_eq; i++)
1688 yif[i] = y0[i];
1689 return (DIFFEQ_ZERO_FOUND);
1690 }
1691
1692 /* adjust step size to stay within interval */
1693 if ((xdiff = xf - *x0) < h_step)
1694 h_step = xdiff;
1695 /* take a step */
1696 x1 = *x0;
1697 rk4_step(y1, x1, y0, dydx0, h_step, n_eq, derivs);
1698 if (stochastic)
1699 /* add stochastic processes */
1700 (*stochastic)(y1, x1, h_step);
1701 x1 += h_step;
1702 /* calculate derivatives and exit function at new point */
1703 (*derivs)(dydx1, y1, x1);
1704 ex1 = (*exit_func)(dydx1, y1, x1);
1705 if (SIGN(ex0) != SIGN(ex1))
1706 break;
1707 /* check for end of interval */
1708 if (FABS(xdiff = xf - x1) < x_accuracy) {
1709 /* end of the interval */
1710 for (i = 0; i < n_eq; i++)
1711 yif[i] = y1[i];
1712 *x0 = x1;
1713 return (DIFFEQ_END_OF_INTERVAL);
1714 }
1715 /* copy the new solution into the old variables */
1716 SWAP_PTR(dydx0, dydx1);
1717 SWAP_PTR(y0, y1);
1718 ex0 = ex1;
1719 *x0 = x1;
1720 } while (1);
1721
1722 if (!exit_func) {
1723 printf("failure in rk_odeint3_na(): solution stepped outside interval\n");
1724 return (DIFFEQ_OUTSIDE_INTERVAL);
1725 }
1726
1727 if (FABS(ex1) < exit_accuracy) {
1728 for (i = 0; i < n_eq; i++)
1729 yif[i] = y1[i];
1730 *x0 = x1;
1731 return (DIFFEQ_ZERO_FOUND);
1732 }
1733
1734 /* The root has been bracketed. */
1735 n_exit_iterations = MAX_EXIT_ITERATIONS;
1736 do {
1737 /* try to take a step to the position where the zero is expected */
1738 /* no stochastic effects are included here! */
1739 h_step = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
1740 x2 = *x0;
1741 rk4_step(y2, x2, y0, dydx0, h_step, n_eq, derivs);
1742 x2 += h_step;
1743 /* check the exit function at the new position */
1744 (*derivs)(dydx2, y2, x2);
1745 ex2 = (*exit_func)(dydx2, y2, x2);
1746 if (FABS(ex2) < exit_accuracy) {
1747 for (i = 0; i < n_eq; i++)
1748 yif[i] = y2[i];
1749 *x0 = x2;
1750 return (DIFFEQ_ZERO_FOUND);
1751 }
1752 /* rebracket the root */
1753 if (SIGN(ex1) == SIGN(ex2)) {
1754 SWAP_PTR(y1, y2);
1755 SWAP_PTR(dydx1, dydx2);
1756 x1 = x2;
1757 ex1 = ex2;
1758 } else {
1759 SWAP_PTR(y0, y2);
1760 SWAP_PTR(dydx0, dydx2);
1761 *x0 = x2;
1762 ex0 = ex2;
1763 }
1764 } while (n_exit_iterations--);
1765 return (DIFFEQ_EXIT_COND_FAILED);
1766}
◆ rk_odeint4()
| long rk_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) ) |
Integrate ODEs until a specific component reaches a target value with intermediate storage.
Integrates a system of ordinary differential equations until a specified component of the solution reaches a target value within a given accuracy. Allows for storing intermediate data points during the integration process.
- Parameters
-
y0 Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Desired accuracies for each variable. accmode Accuracy control modes for each variable. tiny Small values relative to what's important for each variable. misses Array tracking the number of step size reductions per variable. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Desired accuracy for the final value of the independent variable. h_start Suggested starting step size. h_max Maximum step size allowed. h_rec Pointer to store the recommended step size for continuation. exit_value Target value for the specified component. i_exit_value Index of the component to monitor. exit_accuracy Accuracy required for the target value. n_to_skip Number of target value crossings to skip before stopping. store_data Function to store intermediate data points.
- Returns
- DIFFEQ_ZERO_FOUND (>=1) on success, or an error code (<=0) on failure.
Definition at line 1374 of file rkODE.c.
1398 {
1399 double *y_return, *accur;
1400 double *dydx0, *y1, *dydx1, *dydx2, *y2;
1401 double ex0, ex1, ex2, x1, x2, *yscale;
1402 double h_used, h_next, xdiff;
1403 long i, n_exit_iterations, n_step_ups = 0, is_zero;
1404#define MAX_N_STEP_UPS 10
1405
1406 if (*x0 > xf)
1407 return (DIFFEQ_XI_GT_XF);
1408 if (fabs(*x0 - xf) < x_accuracy)
1409 return (DIFFEQ_SOLVED_ALREADY);
1410 if (i_exit_value < 0 || i_exit_value >= n_eq)
1412
1413 /* Meaning of accmode:
1414 * accmode = -1 -> no accuracy control
1415 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
1416 * each step for ith variable. tiny[i] is lower limit
1417 * of significance for the ith variable.
1418 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
1419 * satisfied globally, not locally.
1420 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
1421 * step for the ith variable. tiny[i] is ignored.
1422 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
1423 * be satisfied globally, not locally.
1424 */
1425 for (i = 0; i < n_eq; i++) {
1426 if (accmode[i] < -1 || accmode[i] > 3)
1428 if (accmode[i] < 2 && tiny[i] < TINY)
1429 tiny[i] = TINY;
1430 misses[i] = 0;
1431 }
1432
1433 y_return = y0;
1440
1441 /* calculate derivatives and exit function at the initial point */
1442 (*derivs)(dydx0, y0, *x0);
1443
1444 /* set the scales for evaluating accuracy. yscale[i] is the
1445 * absolute level of accuracy required of the next integration step
1446 */
1448 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1449 accuracy, accur, *x0, xf, n_eq);
1450
1451 ex0 = exit_value - y0[i_exit_value];
1452 if (store_data)
1453 (*store_data)(dydx0, y0, *x0, ex0);
1454 is_zero = 0;
1455 do {
1456 /* check for zero of exit function */
1457 if (fabs(ex0) < exit_accuracy) {
1458 if (!is_zero) {
1459 if (n_to_skip == 0) {
1460 if (store_data)
1461 (*store_data)(dydx0, y0, *x0, ex0);
1462 for (i = 0; i < n_eq; i++)
1463 y_return[i] = y0[i];
1464 *h_rec = h_start;
1465 tfree(dydx0);
1466 tfree(dydx1);
1467 tfree(dydx2);
1468 tfree(yscale);
1469 tfree(accur);
1470 if (y0 != y_return)
1471 tfree(y0);
1472 if (y1 != y_return)
1473 tfree(y1);
1474 if (y2 != y_return)
1475 tfree(y2);
1476 return (DIFFEQ_ZERO_FOUND);
1477 } else {
1478 is_zero = 1;
1479 --n_to_skip;
1480 }
1481 }
1482 } else
1483 is_zero = 0;
1484 /* adjust step size to stay within interval */
1485 if ((xdiff = xf - *x0) < h_start)
1486 h_start = xdiff;
1487 /* take a step */
1488 x1 = *x0;
1489 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
1490 yscale, n_eq, derivs, misses)) {
1491 if (n_step_ups++ > MAX_N_STEP_UPS) {
1493 }
1494 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
1495 continue;
1496 }
1497 /* calculate derivatives and exit function at new point */
1498 (*derivs)(dydx1, y1, x1);
1499 ex1 = exit_value - y1[i_exit_value];
1500 if (store_data)
1501 (*store_data)(dydx1, y1, x1, ex1);
1502 /* check for change in sign of exit function */
1503 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
1504 if (n_to_skip == 0)
1505 break;
1506 else {
1507 --n_to_skip;
1508 is_zero = 1;
1509 }
1510 }
1511 /* check for end of interval */
1512 if (fabs(xdiff = xf - x1) < x_accuracy) {
1513 /* end of the interval */
1514 if (store_data) {
1515 (*derivs)(dydx1, y1, x1);
1516 ex1 = exit_value - y0[i_exit_value];
1517 (*store_data)(dydx1, y1, x1, ex1);
1518 }
1519 for (i = 0; i < n_eq; i++)
1520 y_return[i] = y1[i];
1521 *x0 = x1;
1522 *h_rec = h_start;
1523 tfree(dydx0);
1524 tfree(dydx1);
1525 tfree(dydx2);
1526 tfree(yscale);
1527 tfree(accur);
1528 if (y0 != y_return)
1529 tfree(y0);
1530 if (y1 != y_return)
1531 tfree(y1);
1532 if (y2 != y_return)
1533 tfree(y2);
1534 return (DIFFEQ_END_OF_INTERVAL);
1535 }
1536 /* copy the new solution into the old variables */
1537 SWAP_PTR(dydx0, dydx1);
1538 SWAP_PTR(y0, y1);
1539 ex0 = ex1;
1540 *x0 = x1;
1541 /* adjust the step size as recommended by rk_qcstep() */
1542 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
1543 /* calculate new scale factors */
1544 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1545 accur, n_eq);
1546 } while (1);
1547 *h_rec = h_start;
1548
1549 /* The root has been bracketed. */
1550 n_exit_iterations = MAX_EXIT_ITERATIONS;
1551 do {
1552 /* try to take a step to the position where the zero is expected */
1553 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
1554 x2 = *x0;
1555 /* calculate new scale factors */
1556 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1557 accur, n_eq);
1558 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1559 yscale, n_eq, derivs, misses))
1561 /* check the exit function at the new position */
1562 (*derivs)(dydx2, y2, x2);
1563 ex2 = exit_value - y2[i_exit_value];
1564 if (fabs(ex2) < exit_accuracy) {
1565 for (i = 0; i < n_eq; i++)
1566 y_return[i] = y2[i];
1567 *x0 = x2;
1568 tfree(dydx0);
1569 tfree(dydx1);
1570 tfree(dydx2);
1571 tfree(yscale);
1572 tfree(accur);
1573 if (y0 != y_return)
1574 tfree(y0);
1575 if (y1 != y_return)
1576 tfree(y1);
1577 if (y2 != y_return)
1578 tfree(y2);
1579 return (DIFFEQ_ZERO_FOUND);
1580 }
1581 /* rebracket the root */
1582 if (SIGN(ex1) == SIGN(ex2)) {
1583 SWAP_PTR(y1, y2);
1584 SWAP_PTR(dydx1, dydx2);
1585 x1 = x2;
1586 ex1 = ex2;
1587 } else {
1588 SWAP_PTR(y0, y2);
1589 SWAP_PTR(dydx0, dydx2);
1590 *x0 = x2;
1591 ex0 = ex2;
1592 }
1593 } while (n_exit_iterations--);
1594 return (DIFFEQ_EXIT_COND_FAILED);
1595}
◆ rk_odeint_na()
| long rk_odeint_na | ( | 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, | ||
| double | h_max, | ||
| double * | h_rec ) |
Integrate ODEs without adaptive step-size control.
Integrates a system of ordinary differential equations without using adaptive step-size control. This function is intended for scenarios where adaptive control is not required and allows for faster computations.
- Parameters
-
y0 Initial and final values of dependent variables. derivs Function to compute derivatives. n_eq Number of equations. accuracy Ignored in this function. accmode Ignored in this function. tiny Ignored in this function. misses Ignored in this function. x0 Pointer to the initial value of the independent variable. Updated to final value. xf Upper limit of the independent variable. x_accuracy Ignored in this function. h Step size. h_max Ignored in this function. h_rec Ignored in this function.
- Returns
- DIFFEQ_END_OF_INTERVAL (1) on success, or an error code (<=0) on failure.
Definition at line 1058 of file rkODE.c.
1075 {
1076 double *dydx2, *dydx1, *ytemp, *dydx0;
1077 double xh, hh, h6, x;
1078 long i, j, n;
1079
1080 if ((x = *x0) > xf)
1081 return (DIFFEQ_XI_GT_XF);
1082 if (h == 0)
1083 return (DIFFEQ_ZERO_STEPSIZE);
1084 if (!(n = (xf - x) / h + 0.5))
1085 n = 1;
1086 h = (xf - x) / n;
1087
1092
1093 n_eq--; /* for convenience in loops */
1094
1095 h6 = h / 6;
1096 hh = h / 2;
1097 for (j = 0; j < n; j++) {
1098 if (j == n - 1) {
1099 h = xf - x;
1100 h6 = h / 6;
1101 hh = h / 2;
1102 }
1103 xh = x + hh;
1104
1105 /* first step */
1106 (*derivs)(dydx0, y0, x);
1107 for (i = n_eq; i >= 0; i--)
1108 ytemp[i] = y0[i] + hh * dydx0[i];
1109
1110 /* second step */
1111 (*derivs)(dydx1, ytemp, xh);
1112 for (i = n_eq; i >= 0; i--)
1113 ytemp[i] = y0[i] + hh * dydx1[i];
1114
1115 /* third step */
1116 (*derivs)(dydx2, ytemp, xh);
1117 for (i = n_eq; i >= 0; i--) {
1118 ytemp[i] = y0[i] + h * dydx2[i];
1119 dydx2[i] += dydx1[i];
1120 }
1121
1122 /* fourth step */
1123 (*derivs)(dydx1, ytemp, x = xh + hh);
1124 for (i = n_eq; i >= 0; i--)
1125 y0[i] += h6 * (dydx0[i] + dydx1[i] + 2 * dydx2[i]);
1126 }
1127
1128 tfree(dydx0);
1129 tfree(dydx2);
1130 tfree(dydx1);
1131 tfree(ytemp);
1132 *x0 = x;
1133 return (DIFFEQ_END_OF_INTERVAL);
1134}
◆ rk_qcstep()
| long rk_qcstep | ( | double * | yFinal, |
| double * | x, | ||
| double * | yInitial, | ||
| double * | dydxInitial, | ||
| double | hInput, | ||
| double * | hUsed, | ||
| double * | hRecommended, | ||
| double * | yScale, | ||
| long | equations, | ||
| void(* | derivs )(double *dydx, double *y, double x), | ||
| long * | misses ) |
Definition at line 236 of file rkODE.c.
251 {
252 static MDB_THREAD_LOCAL long last_equations = 0;
253 static MDB_THREAD_LOCAL double *dydxTemp = NULL, *yTemp = NULL;
254 double hOver2, h, xTemp;
255 double error, maxError, hFactor;
256 double localSafetyMargin, localIncreasePower, localDecreasePower, localMaxIncreaseFactor;
257 long i, iWorst = 0, minStepped = 0, noAdaptation = 0;
258
259 mdb_thread_lock(&rk_qctune_lock);
260 localSafetyMargin = safetyMargin;
261 localIncreasePower = increasePower;
262 localDecreasePower = decreasePower;
263 localMaxIncreaseFactor = maxIncreaseFactor;
264 mdb_thread_unlock(&rk_qctune_lock);
265
266 /* for speed, I avoid reallocating the arrays unless it is necessary */
267 if (last_equations < equations) {
268 if (last_equations != 0) {
269 tfree(dydxTemp);
270 tfree(yTemp);
271 }
272 last_equations = equations;
275 }
276
277 for (i = 0; i < equations; i++)
278 if (yScale[i] != DBL_MAX)
279 break;
280 if (i == equations)
281 noAdaptation = 1;
282
283 h = hInput;
284 do {
285#if DEBUG
286 printf("x = %e, h = %e\n", *x, h);
287#endif
288 hOver2 = h / 2;
289 xTemp = *x + hOver2;
290 if (xTemp == *x) {
291 if (!minStepped) {
292 puts("warning: step-size underflow in rk_qcstep()");
293 report_state_dp(stdout, yInitial, dydxInitial, yScale, misses, *x, h, equations);
294 if ((h = 2 * fabs(*x) * DBL_EPSILON) == 0)
295 h = 2 * DBL_EPSILON;
296 hOver2 = h / 2;
297 xTemp = *x + hOver2;
298 minStepped = 1;
299 } else
300 return 0;
301 } else
302 minStepped = 1;
303
304 /* get solution via two half-steps */
305 /* -- first get solution at x+h/2 into yTemp */
306 rk4_step(yTemp, xTemp, yInitial, dydxInitial, hOver2, equations, derivs);
307
308 /* -- get solution at x+h into yFinal, starting from x+h/2 */
309 (*derivs)(dydxTemp, yTemp, xTemp);
310 xTemp = *x + h;
311 rk4_step(yFinal, xTemp, yTemp, dydxTemp, hOver2, equations, derivs);
312
313 /* get solution at x+h into yTemp by taking one step */
314 rk4_step(yTemp, xTemp, yInitial, dydxInitial, h, equations, derivs);
315
316 maxError = 0;
317 for (i = 0; i < equations; i++)
318 yTemp[i] = yFinal[i] - yTemp[i];
319 if (!noAdaptation) {
320 /* Evaluate accuracy */
321 iWorst = -1;
322 for (i = 0; i < equations; i++) {
323#if DEBUG
324 printf("%ld: yTemp=%e yFinal=%e yScale=%e, ", i, yTemp[i], yFinal[i], yScale[i]);
325#endif
326 if (maxError < (error = fabs(yTemp[i]) / yScale[i])) {
327 maxError = error;
328 iWorst = i;
329 }
330#if DEBUG
331 printf(" error = %e\n", error);
332#endif
333 }
334 }
335
336 if (maxError <= 1) {
337 /* step was successful at the following stepsize: */
338 *hUsed = h;
339
340 if (!noAdaptation) {
341 /* Compute recommended stepsize for next step, but don't increase
342 by more than maxIncreaseFactor-fold. Also, use only safetyMargin times
343 the theoretical optimum factor.
344 */
345 if (maxError)
346 hFactor = localSafetyMargin * pow(maxError, -localIncreasePower);
347 else
348 /* error is zero, so go for broke */
349 hFactor = localMaxIncreaseFactor;
350#if DEBUG
351 printf("maxError = %e, hFactor = %e\n", maxError, hFactor);
352#endif
353 if (hFactor > localMaxIncreaseFactor)
354 hFactor = localMaxIncreaseFactor;
355 else if (hFactor < 1)
356 hFactor = 1;
357 } else
358 hFactor = 1;
359
360 *hRecommended = hFactor * h;
361 /* yTemp stores the two-step solution minus the one-step solution.
362 An improved value is obtained by adding 1/15 this difference to
363 the one-step solution .
364 */
365 for (i = 0; i < equations; i++)
366 yFinal[i] += yTemp[i] / 15;
367
368 *x = xTemp; /* update value of independent variable */
369 return (1);
370 }
371
372 if (iWorst >= 0)
373 /* record that equation iWorst caused a step-size reduction */
374 misses[iWorst] += 1;
375
376 /* compute a new, smaller step-size. It will be tested for underflow at the
377 top of the loop.
378 */
379 hOver2 = (h = localSafetyMargin * h * pow(maxError, -localDecreasePower)) / 2;
380 } while (1);
381}
Generated by
1.12.0