Fourth-order Runge-Kutta ODE integration routines (double-precision version). More...
#include "mdb.h"#include <float.h>Go to the source code of this file.
Macros | |
| #define | MAX_EXIT_ITERATIONS 400 |
| #define | ITER_FACTOR 0.995 |
| #define | TINY 1.0e-30 |
| #define | DEBUG 0 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
| #define | MAX_N_STEP_UPS 10 |
Functions | |
| void | new_scale_factors_dp (double *yscale, double *y0, double *dydx0, double h_start, double *tiny, long *accmode, double *accuracy, long n_eq) |
| void | initial_scale_factors_dp (double *yscale, double *y0, double *dydx0, double h_start, double *tiny, long *accmode, double *accuracy, double *accur, double x0, double xf, long n_eq) |
| void | 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. | |
Variables | |
| static double | safetyMargin = 0.9 |
| static double | increasePower = 0.2 |
| static double | decreasePower = 0.25 |
| static double | maxIncreaseFactor = 4.0 |
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.
Macro Definition Documentation
◆ DEBUG
◆ ITER_FACTOR
◆ MAX_EXIT_ITERATIONS
◆ TINY
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 70 of file rkODE.c.
80 {
81 int i;
82
83 for (i = 0; i < n_eq; i++) {
84 if ((accur[i] = accuracy[i]) <= 0) {
85 printf("error: accuracy[%d] = %e (initial_scale_factors_dp)\n", i, accuracy[i]);
86 abort();
87 }
88 switch (accmode[i]) {
89 case -1: /* no accuracy control */
90 yscale[i] = DBL_MAX;
91 break;
92 case 0: /* fractional local accuracy specified */
93 yscale[i] = (y0[i] + dydx0[i] * h_start + tiny[i]) * accur[i];
94 break;
95 case 1: /* fractional global accuracy specified */
96 yscale[i] = (dydx0[i] * h_start + tiny[i]) * accur[i];
97 break;
98 case 2: /* absolute local accuracy specified */
99 yscale[i] = accur[i];
100 break;
101 case 3: /* absolute global accuracy specified */
102 yscale[i] = (accur[i] /= (xf - x0)) * h_start;
103 break;
104 default:
105 printf("error: accmode[%d] = %ld (initial_scale_factors_dp)\n", i, accmode[i]);
106 abort();
107 break;
108 }
109 if (yscale[i] <= 0) {
110 printf("error: yscale[%d] = %e (initial_scale_factors_dp)\n", i, yscale[i]);
111 abort();
112 }
113 }
114}
◆ new_scale_factors_dp()
| void new_scale_factors_dp | ( | double * | yscale, |
| double * | y0, | ||
| double * | dydx0, | ||
| double | h_start, | ||
| double * | tiny, | ||
| long * | accmode, | ||
| double * | accuracy, | ||
| long | n_eq ) |
Definition at line 28 of file rkODE.c.
36 {
37 int i;
38
39 for (i = 0; i < n_eq; i++) {
40 switch (accmode[i]) {
41 case -1: /* no accuracy control */
42 yscale[i] = DBL_MAX;
43 break;
44 case 0: /* fractional local accuracy specified */
45 yscale[i] =
46 (y0[i] + dydx0[i] * h_start + tiny[i]) * accuracy[i];
47 break;
48 case 1: /* fractional global accuracy specified */
49 yscale[i] =
50 (dydx0[i] * h_start + tiny[i]) * accuracy[i];
51 break;
52 case 2: /* absolute local accuracy specified */
53 yscale[i] = accuracy[i];
54 break;
55 case 3: /* absolute global accuracy specified */
56 yscale[i] = accuracy[i] * h_start;
57 break;
58 default:
59 printf("error: accmode[%d] = %ld (new_scale_factors_dp)\n", i, accmode[i]);
60 abort();
61 break;
62 }
63 if (yscale[i] <= 0) {
64 printf("error: yscale[%d] = %e\n", i, yscale[i]);
65 abort();
66 }
67 }
68}
◆ 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 116 of file rkODE.c.
117 {
118 int i;
119 fputs("integration state:\n", fp);
120 fprintf(fp, "%ld equations, indep.var.=%e, step size=%e",
121 n_eq, x, h);
122 fprintf(fp, "\ny : ");
123 for (i = 0; i < n_eq; i++)
124 fprintf(fp, "%e ", y[i]);
125 fprintf(fp, "\ndydx : ");
126 for (i = 0; i < n_eq; i++)
127 fprintf(fp, "%e ", dydx[i]);
128 fprintf(fp, "\ntol.scale: ");
129 for (i = 0; i < n_eq; i++)
130 fprintf(fp, "%e ", yscale[i]);
131 fprintf(fp, "\nmisses : ");
132 for (i = 0; i < n_eq; i++)
133 fprintf(fp, "%ld ", misses[i]);
134}
◆ rk4_qctune()
| void rk4_qctune | ( | double | newSafetyMargin, |
| double | newIncreasePower, | ||
| double | newDecreasePower, | ||
| double | newMaxIncreaseFactor ) |
Definition at line 220 of file rkODE.c.
221 {
222 if (newSafetyMargin > 0 && newSafetyMargin < 1)
223 safetyMargin = newSafetyMargin;
224 if (newIncreasePower > 0)
225 increasePower = newIncreasePower;
226 if (newDecreasePower > 0)
227 decreasePower = newDecreasePower;
228 if (newMaxIncreaseFactor > 1)
229 maxIncreaseFactor = newMaxIncreaseFactor;
230}
◆ 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 141 of file rkODE.c.
150 {
151 static long last_n_eq = 0;
152 static double *k1, *k2, *k3, *yTemp, *dydxTemp;
153 double x1;
154 long i;
155
156 /* for speed, I avoid reallocating the arrays unless it is necessary */
157 if (last_n_eq < n_eq) {
158 if (last_n_eq != 0) {
159 free(k1);
160 free(k2);
161 free(k3);
162 free(yTemp);
163 free(dydx);
164 }
165 last_n_eq = n_eq;
171 }
172
173 /* yf = yi + k1/6 + k2/3 + k3/3 + k4/6 where
174 k1 = h*dydt(x , y)
175 k2 = h*dydt(x+h/2, y+k1/2)
176 k3 = h*dydt(x+h/2, y+k2/2)
177 k4 = h*dydt(x+h , y+k3)
178 */
179
180 /* compute k1 and yTemp for k2 computation */
181 for (i = 0; i < n_eq; i++) {
182 k1[i] = h * dydx[i];
183 yTemp[i] = yi[i] + k1[i] / 2;
184 }
185
186 /* compute k2 and yTemp for k3 computation */
187 x1 = x + h / 2;
188 (*derivs)(dydxTemp, yTemp, x1);
189 for (i = 0; i < n_eq; i++) {
190 k2[i] = h * dydxTemp[i];
191 yTemp[i] = yi[i] + k2[i] / 2;
192 }
193
194 /* compute k3 and yTemp for k4 computation */
195 (*derivs)(dydxTemp, yTemp, x1);
196 for (i = 0; i < n_eq; i++) {
197 k3[i] = h * dydxTemp[i];
198 yTemp[i] = yi[i] + k3[i];
199 }
200
201 /* compute k4 and put results into yf array */
202 x1 = x + h;
203 (*derivs)(dydxTemp, yTemp, x1);
204 for (i = 0; i < n_eq; i++)
205 yf[i] = yi[i] + (k1[i] / 2 + k2[i] + k3[i] + h * dydxTemp[i] / 2) / 3;
206}
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
◆ 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 397 of file rkODE.c.
421 {
422 double *y_return;
423 double *dydx0, *y1, *dydx1, *dydx2, *y2;
424 double ex0, ex1, ex2, *yscale, *accur;
425 double h_used, h_next, x1, x2, xdiff;
426 long i, n_exit_iterations, n_step_ups = 0, is_zero;
427#define MAX_N_STEP_UPS 10
428
429 if (*x0 > xf)
430 return (DIFFEQ_XI_GT_XF);
431 if (FABS(*x0 - xf) < x_accuracy)
432 return (DIFFEQ_SOLVED_ALREADY);
433
434 /* Meaning of accmode:
435 * accmode = -1 -> no accuracy control
436 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
437 * each step for ith variable. tiny[i] is lower limit
438 * of significance for the ith variable.
439 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
440 * satisfied globally, not locally.
441 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
442 * step for the ith variable. tiny[i] is ignored.
443 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
444 * be satisfied globally, not locally.
445 */
446 for (i = 0; i < n_eq; i++) {
447 if (accmode[i] < -1 || accmode[i] > 3)
449 if (accmode[i] < 2 && tiny[i] < TINY)
450 tiny[i] = TINY;
451 misses[i] = 0;
452 }
453
454 y_return = y0;
461
462 /* calculate derivatives and exit function at the initial point */
463 (*derivs)(dydx0, y0, *x0);
464
465 /* set the scales for evaluating accuracy. yscale[i] is the
466 * absolute level of accuracy required of the next integration step.
467 */
469 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
470 accuracy, accur, *x0, xf, n_eq);
471
472 ex0 = exit_func ? (*exit_func)(dydx0, y0, *x0) : 0;
473 if (store_data)
474 (*store_data)(dydx0, y0, *x0, ex0);
475 is_zero = 0;
476
477 do {
478 /* check for zero of exit function */
479 if (exit_func && FABS(ex0) < exit_accuracy) {
480 if (!is_zero) {
481 if (n_to_skip == 0) {
482 if (store_data)
483 (*store_data)(dydx0, y0, *x0, ex0);
484 for (i = 0; i < n_eq; i++)
485 y_return[i] = y0[i];
486 *h_rec = h_start;
487 tfree(dydx0);
488 tfree(dydx1);
489 tfree(dydx2);
490 tfree(yscale);
491 tfree(accur);
492 if (y0 != y_return)
493 tfree(y0);
494 if (y1 != y_return)
495 tfree(y1);
496 if (y2 != y_return)
497 tfree(y2);
498 return (DIFFEQ_ZERO_FOUND);
499 } else {
500 is_zero = 1;
501 --n_to_skip;
502 }
503 }
504 } else
505 is_zero = 0;
506 /* adjust step size to stay within interval */
507 if ((xdiff = xf - *x0) < h_start)
508 h_start = xdiff;
509 /* take a step */
510 x1 = *x0;
511 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
512 yscale, n_eq, derivs, misses)) {
513 if (n_step_ups++ > MAX_N_STEP_UPS)
515 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
516 continue;
517 }
518 /* calculate derivatives and exit function at new point */
519 (*derivs)(dydx1, y1, x1);
520 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
521 if (store_data)
522 (*store_data)(dydx1, y1, x1, ex1);
523 /* check for change in sign of exit function */
524 if (exit_func && SIGN(ex0) != SIGN(ex1) && !is_zero) {
525 if (n_to_skip == 0)
526 break;
527 else {
528 --n_to_skip;
529 is_zero = 1;
530 }
531 }
532
533 /* check for end of interval */
534 if (FABS(xdiff = xf - x1) < x_accuracy) {
535 /* end of the interval */
536 if (store_data) {
537 (*derivs)(dydx1, y1, x1);
538 ex1 = exit_func ? (*exit_func)(dydx1, y1, x1) : 0;
539 (*store_data)(dydx1, y1, x1, ex1);
540 }
541 for (i = 0; i < n_eq; i++)
542 y_return[i] = y1[i];
543 *x0 = x1;
544 *h_rec = h_start;
545 tfree(dydx0);
546 tfree(dydx1);
547 tfree(dydx2);
548 tfree(yscale);
549 tfree(accur);
550 if (y0 != y_return)
551 tfree(y0);
552 if (y1 != y_return)
553 tfree(y1);
554 if (y2 != y_return)
555 tfree(y2);
556 return (DIFFEQ_END_OF_INTERVAL);
557 }
558
559 /* copy the new solution into the old variables */
560 SWAP_PTR(dydx0, dydx1);
561 SWAP_PTR(y0, y1);
562 ex0 = ex1;
563 *x0 = x1;
564 /* adjust the step size as recommended by rk_qcstep() */
565 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
566 /* calculate new scale factors */
567 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
568 accur, n_eq);
569
570 } while (1);
571 *h_rec = h_start;
572
573 if (!exit_func) {
574 printf("failure in rk_odeint(): solution stepped outside interval\n");
575 tfree(dydx0);
576 tfree(dydx1);
577 tfree(dydx2);
578 tfree(yscale);
579 tfree(accur);
580 if (y0 != y_return)
581 tfree(y0);
582 if (y1 != y_return)
583 tfree(y1);
584 if (y2 != y_return)
585 tfree(y2);
586 return (DIFFEQ_OUTSIDE_INTERVAL);
587 }
588
589 /* The root has been bracketed. */
590 n_exit_iterations = MAX_EXIT_ITERATIONS;
591 do {
592 /* try to take a step to the position where the zero is expected */
593 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
594 x2 = *x0;
595 /* calculate new scale factors */
596 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
597 accur, n_eq);
598 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
599 yscale, n_eq, derivs, misses))
601 /* check the exit function at the new position */
602 (*derivs)(dydx2, y2, x2);
603 ex2 = (*exit_func)(dydx2, y2, x2);
604 if (FABS(ex2) < exit_accuracy) {
605 for (i = 0; i < n_eq; i++)
606 y_return[i] = y2[i];
607 *x0 = x2;
608 tfree(dydx0);
609 tfree(dydx1);
610 tfree(dydx2);
611 tfree(yscale);
612 tfree(accur);
613 if (y0 != y_return)
614 tfree(y0);
615 if (y1 != y_return)
616 tfree(y1);
617 if (y2 != y_return)
618 tfree(y2);
619 return (DIFFEQ_ZERO_FOUND);
620 }
621 /* rebracket the root */
622 if (SIGN(ex1) == SIGN(ex2)) {
623 SWAP_PTR(y1, y2);
624 SWAP_PTR(dydx1, dydx2);
625 x1 = x2;
626 ex1 = ex2;
627 } else {
628 SWAP_PTR(y0, y2);
629 SWAP_PTR(dydx0, dydx2);
630 *x0 = x2;
631 ex0 = ex2;
632 }
633 } while (n_exit_iterations--);
634 return (DIFFEQ_EXIT_COND_FAILED);
635}
int tfree(void *ptr)
Frees a memory block and records the deallocation if tracking is enabled.
Definition array.c:230
◆ 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 659 of file rkODE.c.
677 {
678 double *y_return;
679 double *dydx0, *y1, *dydx1, *dydx2, *y2;
680 double *yscale, *accur;
681 double x1;
682 double h_used, h_next, xdiff;
683 long i, n_step_ups = 0;
684#define MAX_N_STEP_UPS 10
685
686 if (*x0 > xf)
687 return (DIFFEQ_XI_GT_XF);
688 if (FABS(*x0 - xf) < x_accuracy)
689 return (DIFFEQ_SOLVED_ALREADY);
690
691 /* Meaning of accmode:
692 * accmode = -1 -> no accuracy control
693 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
694 * each step for ith variable. tiny[i] is lower limit
695 * of significance for the ith variable.
696 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
697 * satisfied globally, not locally.
698 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
699 * step for the ith variable. tiny[i] is ignored.
700 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
701 * be satisfied globally, not locally.
702 */
703 for (i = 0; i < n_eq; i++) {
704 if (accmode[i] < -1 || accmode[i] > 3)
706 if (accmode[i] < 2 && tiny[i] < TINY)
707 tiny[i] = TINY;
708 misses[i] = 0;
709 }
710
711 y_return = y0;
718
719 /* calculate derivatives at the initial point */
720 (*derivs)(dydx0, y0, *x0);
721
722 /* set the scales for evaluating accuracy. yscale[i] is the
723 * absolute level of accuracy required of the next integration step
724 */
726 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
727 accuracy, accur, *x0, xf, n_eq);
728
729 do {
730 /* adjust step size to stay within interval */
731 if ((xdiff = xf - *x0) < h_start)
732 h_start = xdiff;
733 /* take a step */
734 x1 = *x0;
735 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
736 yscale, n_eq, derivs, misses)) {
737 if (n_step_ups++ > MAX_N_STEP_UPS) {
738 puts("error: cannot take step (rk_odeint1--1)");
739 printf("xf = %.16e x0 = %.16e\n", xf, *x0);
740 printf("h_start = %.16e h_used = %.16e\n", h_start, h_used);
741 puts("dump of integration state:");
742 puts(" variable value derivative scale misses");
743 puts("---------------------------------------------------------------");
744 for (i = 0; i < n_eq; i++)
745 printf(" %5ld %13.6e %13.6e %13.6e %5ld \n",
746 i, y0[i], dydx0[i], yscale[i], misses[i]);
747 exit(1);
748 }
749 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
750 continue;
751 }
752 /* check for end of interval */
753 if (FABS(xdiff = xf - x1) < x_accuracy) {
754 /* end of the interval */
755 for (i = 0; i < n_eq; i++)
756 y_return[i] = y1[i];
757 *x0 = x1;
758 *h_rec = h_start;
759 tfree(dydx0);
760 tfree(dydx1);
761 tfree(dydx2);
762 tfree(yscale);
763 tfree(accur);
764 if (y0 != y_return)
765 tfree(y0);
766 if (y1 != y_return)
767 tfree(y1);
768 if (y2 != y_return)
769 tfree(y2);
770 return (DIFFEQ_END_OF_INTERVAL);
771 }
772 /* calculate derivatives at new point */
773 (*derivs)(dydx1, y1, x1);
774 /* copy the new solution into the old variables */
775 SWAP_PTR(dydx0, dydx1);
776 SWAP_PTR(y0, y1);
777 *x0 = x1;
778 /* adjust the step size as recommended by rk_qcstep() */
779 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
780 /* calculate new scale factors */
781 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
782 accur, n_eq);
783 } while (1);
784}
◆ 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 812 of file rkODE.c.
835 {
836 double *y_return, *accur;
837 double *dydx0, *y1, *dydx1, *dydx2, *y2;
838 double ex0, ex1, ex2, x1, x2, *yscale;
839 double h_used, h_next, xdiff;
840 long i, n_exit_iterations, n_step_ups = 0, is_zero;
841#define MAX_N_STEP_UPS 10
842
843 if (*x0 > xf)
844 return (DIFFEQ_XI_GT_XF);
845 if (FABS(*x0 - xf) < x_accuracy)
846 return (DIFFEQ_SOLVED_ALREADY);
847 if (i_exit_value < 0 || i_exit_value >= n_eq)
849
850 /* Meaning of accmode:
851 * accmode = -1 -> no accuracy control
852 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
853 * each step for ith variable. tiny[i] is lower limit
854 * of significance for the ith variable.
855 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
856 * satisfied globally, not locally.
857 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
858 * step for the ith variable. tiny[i] is ignored.
859 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
860 * be satisfied globally, not locally.
861 */
862 for (i = 0; i < n_eq; i++) {
863 if (accmode[i] < -1 || accmode[i] > 3)
865 if (accmode[i] < 2 && tiny[i] < TINY)
866 tiny[i] = TINY;
867 misses[i] = 0;
868 }
869
870 y_return = y0;
877
878 /* calculate derivatives and exit function at the initial point */
879 (*derivs)(dydx0, y0, *x0);
880
881 /* set the scales for evaluating accuracy. yscale[i] is the
882 * absolute level of accuracy required of the next integration step
883 */
885 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
886 accuracy, accur, *x0, xf, n_eq);
887
888 ex0 = exit_value - y0[i_exit_value];
889 is_zero = 0;
890 do {
891 /* check for zero of exit function */
892 if (FABS(ex0) < exit_accuracy) {
893 if (!is_zero) {
894 if (n_to_skip == 0) {
895 for (i = 0; i < n_eq; i++)
896 y_return[i] = y0[i];
897 *h_rec = h_start;
898 tfree(dydx0);
899 tfree(dydx1);
900 tfree(dydx2);
901 tfree(yscale);
902 tfree(accur);
903 if (y0 != y_return)
904 tfree(y0);
905 if (y1 != y_return)
906 tfree(y1);
907 if (y2 != y_return)
908 tfree(y2);
909 return (DIFFEQ_ZERO_FOUND);
910 } else {
911 is_zero = 1;
912 --n_to_skip;
913 }
914 }
915 } else
916 is_zero = 0;
917 /* adjust step size to stay within interval */
918 if ((xdiff = xf - *x0) < h_start)
919 h_start = xdiff;
920 /* take a step */
921 x1 = *x0;
922 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
923 yscale, n_eq, derivs, misses)) {
924 if (n_step_ups++ > MAX_N_STEP_UPS) {
926 }
927 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
928 continue;
929 }
930 /* calculate derivatives and exit function at new point */
931 (*derivs)(dydx1, y1, x1);
932 ex1 = exit_value - y1[i_exit_value];
933 /* check for change in sign of exit function */
934 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
935 if (n_to_skip == 0)
936 break;
937 else {
938 --n_to_skip;
939 is_zero = 1;
940 }
941 }
942 /* check for end of interval */
943 if (FABS(xdiff = xf - x1) < x_accuracy) {
944 /* end of the interval */
945 for (i = 0; i < n_eq; i++)
946 y_return[i] = y1[i];
947 *x0 = x1;
948 *h_rec = h_start;
949 tfree(dydx0);
950 tfree(dydx1);
951 tfree(dydx2);
952 tfree(yscale);
953 tfree(accur);
954 if (y0 != y_return)
955 tfree(y0);
956 if (y1 != y_return)
957 tfree(y1);
958 if (y2 != y_return)
959 tfree(y2);
960 return (DIFFEQ_END_OF_INTERVAL);
961 }
962 /* copy the new solution into the old variables */
963 SWAP_PTR(dydx0, dydx1);
964 SWAP_PTR(y0, y1);
965 ex0 = ex1;
966 *x0 = x1;
967 /* adjust the step size as recommended by rk_qcstep() */
968 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
969 /* calculate new scale factors */
970 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
971 accur, n_eq);
972
973 } while (1);
974 *h_rec = h_start;
975
976 /* The root has been bracketed. */
977 n_exit_iterations = MAX_EXIT_ITERATIONS;
978 do {
979 /* try to take a step to the position where the zero is expected */
980 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
981 x2 = *x0;
982 /* calculate new scale factors */
983 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
984 accur, n_eq);
985 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
986 yscale, n_eq, derivs, misses))
988 /* check the exit function at the new position */
989 (*derivs)(dydx2, y2, x2);
990 ex2 = exit_value - y2[i_exit_value];
991 if (FABS(ex2) < exit_accuracy) {
992 for (i = 0; i < n_eq; i++)
993 y_return[i] = y2[i];
994 *x0 = x2;
995 tfree(dydx0);
996 tfree(dydx1);
997 tfree(dydx2);
998 tfree(yscale);
999 tfree(accur);
1000 if (y0 != y_return)
1001 tfree(y0);
1002 if (y1 != y_return)
1003 tfree(y1);
1004 if (y2 != y_return)
1005 tfree(y2);
1006 return (DIFFEQ_ZERO_FOUND);
1007 }
1008 /* rebracket the root */
1009 if (SIGN(ex1) == SIGN(ex2)) {
1010 SWAP_PTR(y1, y2);
1011 SWAP_PTR(dydx1, dydx2);
1012 x1 = x2;
1013 ex1 = ex2;
1014 } else {
1015 SWAP_PTR(y0, y2);
1016 SWAP_PTR(dydx0, dydx2);
1017 *x0 = x2;
1018 ex0 = ex2;
1019 }
1020 } while (n_exit_iterations--);
1021 return (DIFFEQ_EXIT_COND_FAILED);
1022}
◆ 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 1147 of file rkODE.c.
1168 {
1169 static double *y0, *yscale;
1170 static double *dydx0, *y1, *dydx1, *dydx2, *y2, *accur;
1171 static long last_neq = 0;
1172 double ex0, ex1, ex2, x1, x2;
1173 double h_used, h_next, xdiff;
1174 long i, n_exit_iterations, n_step_ups = 0;
1175#define MAX_N_STEP_UPS 10
1176
1177 if (*x0 > xf)
1178 return (DIFFEQ_XI_GT_XF);
1179 if (FABS(*x0 - xf) < x_accuracy)
1180 return (DIFFEQ_SOLVED_ALREADY);
1181
1182 /* Meaning of accmode:
1183 * accmode = -1 -> no accuracy control
1184 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
1185 * each step for ith variable. tiny[i] is lower limit
1186 * of significance for the ith variable.
1187 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
1188 * satisfied globally, not locally.
1189 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
1190 * step for the ith variable. tiny[i] is ignored.
1191 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
1192 * be satisfied globally, not locally.
1193 */
1194 for (i = 0; i < n_eq; i++) {
1195 if (accmode[i] < -1 || accmode[i] > 3)
1197 if (accmode[i] < 2 && tiny[i] < TINY)
1198 tiny[i] = TINY;
1199 misses[i] = 0;
1200 }
1201
1202 if (last_neq < n_eq) {
1203 if (last_neq != 0) {
1204 tfree(y0);
1205 tfree(dydx0);
1206 tfree(y1);
1207 tfree(dydx1);
1208 tfree(y2);
1209 tfree(dydx2);
1210 tfree(yscale);
1211 tfree(accur);
1212 }
1221 last_neq = n_eq;
1222 }
1223
1224 for (i = 0; i < n_eq; i++)
1225 y0[i] = yif[i];
1226
1227 /* calculate derivatives and exit function at the initial point */
1228 (*derivs)(dydx0, y0, *x0);
1229
1230 /* set the scales for evaluating accuracy. yscale[i] is the
1231 * absolute level of accuracy required of the next integration step
1232 */
1233 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1234 accuracy, accur, *x0, xf, n_eq);
1235
1236 ex0 = (*exit_func)(dydx0, y0, *x0);
1237
1238 do {
1239 /* check for zero of exit function */
1240 if (FABS(ex0) < exit_accuracy) {
1241 for (i = 0; i < n_eq; i++)
1242 yif[i] = y0[i];
1243 *h_rec = h_start;
1244 return (DIFFEQ_ZERO_FOUND);
1245 }
1246
1247 /* adjust step size to stay within interval */
1248 if ((xdiff = xf - *x0) < h_start)
1249 h_start = xdiff;
1250 /* take a step */
1251 x1 = *x0;
1252 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
1253 yscale, n_eq, derivs, misses)) {
1254 if (n_step_ups++ > MAX_N_STEP_UPS)
1256 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
1257 continue;
1258 }
1259 /* calculate derivatives and exit function at new point */
1260 (*derivs)(dydx1, y1, x1);
1261 ex1 = (*exit_func)(dydx1, y1, x1);
1262 if (SIGN(ex0) != SIGN(ex1))
1263 break;
1264 /* check for end of interval */
1265 if (FABS(xdiff = xf - x1) < x_accuracy) {
1266 /* end of the interval */
1267 for (i = 0; i < n_eq; i++)
1268 yif[i] = y1[i];
1269 *x0 = x1;
1270 *h_rec = h_start;
1271 return (DIFFEQ_END_OF_INTERVAL);
1272 }
1273 /* copy the new solution into the old variables */
1274 SWAP_PTR(dydx0, dydx1);
1275 SWAP_PTR(y0, y1);
1276 ex0 = ex1;
1277 *x0 = x1;
1278 /* adjust the step size as recommended by rk_qcstep() */
1279 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
1280 /* calculate new scale factors */
1281 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1282 accur, n_eq);
1283 } while (1);
1284 *h_rec = h_start;
1285
1286 if (!exit_func) {
1287 printf("failure in rk_odeint3(): solution stepped outside interval\n");
1288 return (DIFFEQ_OUTSIDE_INTERVAL);
1289 }
1290
1291 if (FABS(ex1) < exit_accuracy) {
1292 for (i = 0; i < n_eq; i++)
1293 yif[i] = y1[i];
1294 *x0 = x1;
1295 return (DIFFEQ_ZERO_FOUND);
1296 }
1297
1298 /* The root has been bracketed. */
1299 n_exit_iterations = MAX_EXIT_ITERATIONS;
1300 do {
1301 /* try to take a step to the position where the zero is expected */
1302 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
1303 x2 = *x0;
1304 /* calculate new scale factors */
1305 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1306 accur, n_eq);
1307 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1308 yscale, n_eq, derivs, misses))
1310 /* check the exit function at the new position */
1311 (*derivs)(dydx2, y2, x2);
1312 ex2 = (*exit_func)(dydx2, y2, x2);
1313 if (FABS(ex2) < exit_accuracy) {
1314 for (i = 0; i < n_eq; i++)
1315 yif[i] = y2[i];
1316 *x0 = x2;
1317 return (DIFFEQ_ZERO_FOUND);
1318 }
1319 /* rebracket the root */
1320 if (SIGN(ex1) == SIGN(ex2)) {
1321 SWAP_PTR(y1, y2);
1322 SWAP_PTR(dydx1, dydx2);
1323 x1 = x2;
1324 ex1 = ex2;
1325 } else {
1326 SWAP_PTR(y0, y2);
1327 SWAP_PTR(dydx0, dydx2);
1328 *x0 = x2;
1329 ex0 = ex2;
1330 }
1331 } while (n_exit_iterations--);
1332 return (DIFFEQ_EXIT_COND_FAILED);
1333}
◆ 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 1609 of file rkODE.c.
1630 {
1631 static double *y0, *yscale;
1632 static double *dydx0, *y1, *dydx1, *dydx2, *y2, *accur;
1633 static long last_neq = 0;
1634 double ex0, ex1, ex2, x1, x2;
1635 double xdiff;
1636 long i, n_exit_iterations;
1637#define MAX_N_STEP_UPS 10
1638
1639 if (*x0 > xf)
1640 return (DIFFEQ_XI_GT_XF);
1641 if (FABS(*x0 - xf) < x_accuracy)
1642 return (DIFFEQ_SOLVED_ALREADY);
1643
1644 if (last_neq < n_eq) {
1645 if (last_neq != 0) {
1646 tfree(y0);
1647 tfree(dydx0);
1648 tfree(y1);
1649 tfree(dydx1);
1650 tfree(y2);
1651 tfree(dydx2);
1652 tfree(yscale);
1653 tfree(accur);
1654 }
1661 last_neq = n_eq;
1662 }
1663
1664 for (i = 0; i < n_eq; i++)
1665 y0[i] = yif[i];
1666
1667 /* calculate derivatives and exit function at the initial point */
1668 (*derivs)(dydx0, y0, *x0);
1669
1670 ex0 = (*exit_func)(dydx0, y0, *x0);
1671
1672 do {
1673 /* check for zero of exit function */
1674 if (FABS(ex0) < exit_accuracy) {
1675 for (i = 0; i < n_eq; i++)
1676 yif[i] = y0[i];
1677 return (DIFFEQ_ZERO_FOUND);
1678 }
1679
1680 /* adjust step size to stay within interval */
1681 if ((xdiff = xf - *x0) < h_step)
1682 h_step = xdiff;
1683 /* take a step */
1684 x1 = *x0;
1685 rk4_step(y1, x1, y0, dydx0, h_step, n_eq, derivs);
1686 if (stochastic)
1687 /* add stochastic processes */
1688 (*stochastic)(y1, x1, h_step);
1689 x1 += h_step;
1690 /* calculate derivatives and exit function at new point */
1691 (*derivs)(dydx1, y1, x1);
1692 ex1 = (*exit_func)(dydx1, y1, x1);
1693 if (SIGN(ex0) != SIGN(ex1))
1694 break;
1695 /* check for end of interval */
1696 if (FABS(xdiff = xf - x1) < x_accuracy) {
1697 /* end of the interval */
1698 for (i = 0; i < n_eq; i++)
1699 yif[i] = y1[i];
1700 *x0 = x1;
1701 return (DIFFEQ_END_OF_INTERVAL);
1702 }
1703 /* copy the new solution into the old variables */
1704 SWAP_PTR(dydx0, dydx1);
1705 SWAP_PTR(y0, y1);
1706 ex0 = ex1;
1707 *x0 = x1;
1708 } while (1);
1709
1710 if (!exit_func) {
1711 printf("failure in rk_odeint3_na(): solution stepped outside interval\n");
1712 return (DIFFEQ_OUTSIDE_INTERVAL);
1713 }
1714
1715 if (FABS(ex1) < exit_accuracy) {
1716 for (i = 0; i < n_eq; i++)
1717 yif[i] = y1[i];
1718 *x0 = x1;
1719 return (DIFFEQ_ZERO_FOUND);
1720 }
1721
1722 /* The root has been bracketed. */
1723 n_exit_iterations = MAX_EXIT_ITERATIONS;
1724 do {
1725 /* try to take a step to the position where the zero is expected */
1726 /* no stochastic effects are included here! */
1727 h_step = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
1728 x2 = *x0;
1729 rk4_step(y2, x2, y0, dydx0, h_step, n_eq, derivs);
1730 x2 += h_step;
1731 /* check the exit function at the new position */
1732 (*derivs)(dydx2, y2, x2);
1733 ex2 = (*exit_func)(dydx2, y2, x2);
1734 if (FABS(ex2) < exit_accuracy) {
1735 for (i = 0; i < n_eq; i++)
1736 yif[i] = y2[i];
1737 *x0 = x2;
1738 return (DIFFEQ_ZERO_FOUND);
1739 }
1740 /* rebracket the root */
1741 if (SIGN(ex1) == SIGN(ex2)) {
1742 SWAP_PTR(y1, y2);
1743 SWAP_PTR(dydx1, dydx2);
1744 x1 = x2;
1745 ex1 = ex2;
1746 } else {
1747 SWAP_PTR(y0, y2);
1748 SWAP_PTR(dydx0, dydx2);
1749 *x0 = x2;
1750 ex0 = ex2;
1751 }
1752 } while (n_exit_iterations--);
1753 return (DIFFEQ_EXIT_COND_FAILED);
1754}
◆ 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 1362 of file rkODE.c.
1386 {
1387 double *y_return, *accur;
1388 double *dydx0, *y1, *dydx1, *dydx2, *y2;
1389 double ex0, ex1, ex2, x1, x2, *yscale;
1390 double h_used, h_next, xdiff;
1391 long i, n_exit_iterations, n_step_ups = 0, is_zero;
1392#define MAX_N_STEP_UPS 10
1393
1394 if (*x0 > xf)
1395 return (DIFFEQ_XI_GT_XF);
1396 if (fabs(*x0 - xf) < x_accuracy)
1397 return (DIFFEQ_SOLVED_ALREADY);
1398 if (i_exit_value < 0 || i_exit_value >= n_eq)
1400
1401 /* Meaning of accmode:
1402 * accmode = -1 -> no accuracy control
1403 * accmode = 0 -> accuracy[i] is desired fractional accuracy at
1404 * each step for ith variable. tiny[i] is lower limit
1405 * of significance for the ith variable.
1406 * accmode = 1 -> same as accmode=0, except that the accuracy is to be
1407 * satisfied globally, not locally.
1408 * accmode = 2 -> accuracy[i] is the desired absolute accuracy per
1409 * step for the ith variable. tiny[i] is ignored.
1410 * accmode = 3 -> samed as accmode=2, except that the accuracy is to
1411 * be satisfied globally, not locally.
1412 */
1413 for (i = 0; i < n_eq; i++) {
1414 if (accmode[i] < -1 || accmode[i] > 3)
1416 if (accmode[i] < 2 && tiny[i] < TINY)
1417 tiny[i] = TINY;
1418 misses[i] = 0;
1419 }
1420
1421 y_return = y0;
1428
1429 /* calculate derivatives and exit function at the initial point */
1430 (*derivs)(dydx0, y0, *x0);
1431
1432 /* set the scales for evaluating accuracy. yscale[i] is the
1433 * absolute level of accuracy required of the next integration step
1434 */
1436 initial_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1437 accuracy, accur, *x0, xf, n_eq);
1438
1439 ex0 = exit_value - y0[i_exit_value];
1440 if (store_data)
1441 (*store_data)(dydx0, y0, *x0, ex0);
1442 is_zero = 0;
1443 do {
1444 /* check for zero of exit function */
1445 if (fabs(ex0) < exit_accuracy) {
1446 if (!is_zero) {
1447 if (n_to_skip == 0) {
1448 if (store_data)
1449 (*store_data)(dydx0, y0, *x0, ex0);
1450 for (i = 0; i < n_eq; i++)
1451 y_return[i] = y0[i];
1452 *h_rec = h_start;
1453 tfree(dydx0);
1454 tfree(dydx1);
1455 tfree(dydx2);
1456 tfree(yscale);
1457 tfree(accur);
1458 if (y0 != y_return)
1459 tfree(y0);
1460 if (y1 != y_return)
1461 tfree(y1);
1462 if (y2 != y_return)
1463 tfree(y2);
1464 return (DIFFEQ_ZERO_FOUND);
1465 } else {
1466 is_zero = 1;
1467 --n_to_skip;
1468 }
1469 }
1470 } else
1471 is_zero = 0;
1472 /* adjust step size to stay within interval */
1473 if ((xdiff = xf - *x0) < h_start)
1474 h_start = xdiff;
1475 /* take a step */
1476 x1 = *x0;
1477 if (!rk_qcstep(y1, &x1, y0, dydx0, h_start, &h_used, &h_next,
1478 yscale, n_eq, derivs, misses)) {
1479 if (n_step_ups++ > MAX_N_STEP_UPS) {
1481 }
1482 h_start = (n_step_ups - 1 ? h_start * 10 : h_used * 10);
1483 continue;
1484 }
1485 /* calculate derivatives and exit function at new point */
1486 (*derivs)(dydx1, y1, x1);
1487 ex1 = exit_value - y1[i_exit_value];
1488 if (store_data)
1489 (*store_data)(dydx1, y1, x1, ex1);
1490 /* check for change in sign of exit function */
1491 if (SIGN(ex0) != SIGN(ex1) && !is_zero) {
1492 if (n_to_skip == 0)
1493 break;
1494 else {
1495 --n_to_skip;
1496 is_zero = 1;
1497 }
1498 }
1499 /* check for end of interval */
1500 if (fabs(xdiff = xf - x1) < x_accuracy) {
1501 /* end of the interval */
1502 if (store_data) {
1503 (*derivs)(dydx1, y1, x1);
1504 ex1 = exit_value - y0[i_exit_value];
1505 (*store_data)(dydx1, y1, x1, ex1);
1506 }
1507 for (i = 0; i < n_eq; i++)
1508 y_return[i] = y1[i];
1509 *x0 = x1;
1510 *h_rec = h_start;
1511 tfree(dydx0);
1512 tfree(dydx1);
1513 tfree(dydx2);
1514 tfree(yscale);
1515 tfree(accur);
1516 if (y0 != y_return)
1517 tfree(y0);
1518 if (y1 != y_return)
1519 tfree(y1);
1520 if (y2 != y_return)
1521 tfree(y2);
1522 return (DIFFEQ_END_OF_INTERVAL);
1523 }
1524 /* copy the new solution into the old variables */
1525 SWAP_PTR(dydx0, dydx1);
1526 SWAP_PTR(y0, y1);
1527 ex0 = ex1;
1528 *x0 = x1;
1529 /* adjust the step size as recommended by rk_qcstep() */
1530 h_start = (h_next > h_max ? (h_max ? h_max : h_next) : h_next);
1531 /* calculate new scale factors */
1532 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1533 accur, n_eq);
1534 } while (1);
1535 *h_rec = h_start;
1536
1537 /* The root has been bracketed. */
1538 n_exit_iterations = MAX_EXIT_ITERATIONS;
1539 do {
1540 /* try to take a step to the position where the zero is expected */
1541 h_start = -ex0 * (x1 - *x0) / (ex1 - ex0) * ITER_FACTOR;
1542 x2 = *x0;
1543 /* calculate new scale factors */
1544 new_scale_factors_dp(yscale, y0, dydx0, h_start, tiny, accmode,
1545 accur, n_eq);
1546 if (!rk_qcstep(y2, &x2, y0, dydx0, h_start, &h_used, &h_next,
1547 yscale, n_eq, derivs, misses))
1549 /* check the exit function at the new position */
1550 (*derivs)(dydx2, y2, x2);
1551 ex2 = exit_value - y2[i_exit_value];
1552 if (fabs(ex2) < exit_accuracy) {
1553 for (i = 0; i < n_eq; i++)
1554 y_return[i] = y2[i];
1555 *x0 = x2;
1556 tfree(dydx0);
1557 tfree(dydx1);
1558 tfree(dydx2);
1559 tfree(yscale);
1560 tfree(accur);
1561 if (y0 != y_return)
1562 tfree(y0);
1563 if (y1 != y_return)
1564 tfree(y1);
1565 if (y2 != y_return)
1566 tfree(y2);
1567 return (DIFFEQ_ZERO_FOUND);
1568 }
1569 /* rebracket the root */
1570 if (SIGN(ex1) == SIGN(ex2)) {
1571 SWAP_PTR(y1, y2);
1572 SWAP_PTR(dydx1, dydx2);
1573 x1 = x2;
1574 ex1 = ex2;
1575 } else {
1576 SWAP_PTR(y0, y2);
1577 SWAP_PTR(dydx0, dydx2);
1578 *x0 = x2;
1579 ex0 = ex2;
1580 }
1581 } while (n_exit_iterations--);
1582 return (DIFFEQ_EXIT_COND_FAILED);
1583}
◆ 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 1046 of file rkODE.c.
1063 {
1064 double *dydx2, *dydx1, *ytemp, *dydx0;
1065 double xh, hh, h6, x;
1066 long i, j, n;
1067
1068 if ((x = *x0) > xf)
1069 return (DIFFEQ_XI_GT_XF);
1070 if (h == 0)
1071 return (DIFFEQ_ZERO_STEPSIZE);
1072 if (!(n = (xf - x) / h + 0.5))
1073 n = 1;
1074 h = (xf - x) / n;
1075
1080
1081 n_eq--; /* for convenience in loops */
1082
1083 h6 = h / 6;
1084 hh = h / 2;
1085 for (j = 0; j < n; j++) {
1086 if (j == n - 1) {
1087 h = xf - x;
1088 h6 = h / 6;
1089 hh = h / 2;
1090 }
1091 xh = x + hh;
1092
1093 /* first step */
1094 (*derivs)(dydx0, y0, x);
1095 for (i = n_eq; i >= 0; i--)
1096 ytemp[i] = y0[i] + hh * dydx0[i];
1097
1098 /* second step */
1099 (*derivs)(dydx1, ytemp, xh);
1100 for (i = n_eq; i >= 0; i--)
1101 ytemp[i] = y0[i] + hh * dydx1[i];
1102
1103 /* third step */
1104 (*derivs)(dydx2, ytemp, xh);
1105 for (i = n_eq; i >= 0; i--) {
1106 ytemp[i] = y0[i] + h * dydx2[i];
1107 dydx2[i] += dydx1[i];
1108 }
1109
1110 /* fourth step */
1111 (*derivs)(dydx1, ytemp, x = xh + hh);
1112 for (i = n_eq; i >= 0; i--)
1113 y0[i] += h6 * (dydx0[i] + dydx1[i] + 2 * dydx2[i]);
1114 }
1115
1116 tfree(dydx0);
1117 tfree(dydx2);
1118 tfree(dydx1);
1119 tfree(ytemp);
1120 *x0 = x;
1121 return (DIFFEQ_END_OF_INTERVAL);
1122}
◆ 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 232 of file rkODE.c.
247 {
248 static long last_equations = 0;
249 static double *dydxTemp, *yTemp;
250 double hOver2, h, xTemp;
251 double error, maxError, hFactor;
252 long i, iWorst = 0, minStepped = 0, noAdaptation = 0;
253
254 /* for speed, I avoid reallocating the arrays unless it is necessary */
255 if (last_equations < equations) {
256 if (last_equations != 0) {
257 tfree(dydxTemp);
258 tfree(yTemp);
259 }
260 last_equations = equations;
263 }
264
265 for (i = 0; i < equations; i++)
266 if (yScale[i] != DBL_MAX)
267 break;
268 if (i == equations)
269 noAdaptation = 1;
270
271 h = hInput;
272 do {
273#if DEBUG
274 printf("x = %e, h = %e\n", *x, h);
275#endif
276 hOver2 = h / 2;
277 xTemp = *x + hOver2;
278 if (xTemp == *x) {
279 if (!minStepped) {
280 puts("warning: step-size underflow in rk_qcstep()");
281 report_state_dp(stdout, yInitial, dydxInitial, yScale, misses, *x, h, equations);
282 if ((h = 2 * fabs(*x) * DBL_EPSILON) == 0)
283 h = 2 * DBL_EPSILON;
284 hOver2 = h / 2;
285 xTemp = *x + hOver2;
286 minStepped = 1;
287 } else
288 return 0;
289 } else
290 minStepped = 1;
291
292 /* get solution via two half-steps */
293 /* -- first get solution at x+h/2 into yTemp */
294 rk4_step(yTemp, xTemp, yInitial, dydxInitial, hOver2, equations, derivs);
295
296 /* -- get solution at x+h into yFinal, starting from x+h/2 */
297 (*derivs)(dydxTemp, yTemp, xTemp);
298 xTemp = *x + h;
299 rk4_step(yFinal, xTemp, yTemp, dydxTemp, hOver2, equations, derivs);
300
301 /* get solution at x+h into yTemp by taking one step */
302 rk4_step(yTemp, xTemp, yInitial, dydxInitial, h, equations, derivs);
303
304 maxError = 0;
305 for (i = 0; i < equations; i++)
306 yTemp[i] = yFinal[i] - yTemp[i];
307 if (!noAdaptation) {
308 /* Evaluate accuracy */
309 iWorst = -1;
310 for (i = 0; i < equations; i++) {
311#if DEBUG
312 printf("%ld: yTemp=%e yFinal=%e yScale=%e, ", i, yTemp[i], yFinal[i], yScale[i]);
313#endif
314 if (maxError < (error = fabs(yTemp[i]) / yScale[i])) {
315 maxError = error;
316 iWorst = i;
317 }
318#if DEBUG
319 printf(" error = %e\n", error);
320#endif
321 }
322 }
323
324 if (maxError <= 1) {
325 /* step was successful at the following stepsize: */
326 *hUsed = h;
327
328 if (!noAdaptation) {
329 /* Compute recommended stepsize for next step, but don't increase
330 by more than maxIncreaseFactor-fold. Also, use only safetyMargin times
331 the theoretical optimum factor.
332 */
333 if (maxError)
334 hFactor = safetyMargin * pow(maxError, -increasePower);
335 else
336 /* error is zero, so go for broke */
337 hFactor = maxIncreaseFactor;
338#if DEBUG
339 printf("maxError = %e, hFactor = %e\n", maxError, hFactor);
340#endif
341 if (hFactor > maxIncreaseFactor)
342 hFactor = maxIncreaseFactor;
343 else if (hFactor < 1)
344 hFactor = 1;
345 } else
346 hFactor = 1;
347
348 *hRecommended = hFactor * h;
349 /* yTemp stores the two-step solution minus the one-step solution.
350 An improved value is obtained by adding 1/15 this difference to
351 the one-step solution .
352 */
353 for (i = 0; i < equations; i++)
354 yFinal[i] += yTemp[i] / 15;
355
356 *x = xTemp; /* update value of independent variable */
357 return (1);
358 }
359
360 if (iWorst >= 0)
361 /* record that equation iWorst caused a step-size reduction */
362 misses[iWorst] += 1;
363
364 /* compute a new, smaller step-size. It will be tested for underflow at the
365 top of the loop.
366 */
367 hOver2 = (h = safetyMargin * h * pow(maxError, -decreasePower)) / 2;
368 } while (1);
369}
Variable Documentation
◆ decreasePower
◆ increasePower
◆ maxIncreaseFactor
◆ safetyMargin
Generated by
1.12.0