|
|
◆ dopn43_r()
| void dopn43_r |
( |
int |
n, |
|
|
double * |
t, |
|
|
double |
y[], |
|
|
double |
yp[], |
|
|
double |
tout, |
|
|
double |
tend, |
|
|
double * |
rtol, |
|
|
double * |
atol, |
|
|
int |
itol, |
|
|
int |
mode, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
int |
iwork[], |
|
|
int |
liwork, |
|
|
int * |
info, |
|
|
double * |
tt, |
|
|
double |
yy[], |
|
|
double |
yypp[], |
|
|
int * |
irev |
|
) |
| |
Initial value problem of ordinary differential equations (4(3)-th order Runge-Kutta-Nystrom method) (for second order differential equations) (Reverse communication version)
- Purpose
- This program integrates a system of second order ordinary differential equations of the form
d2y/dt2 = f(t, y), y = y0, y' = y'0 at t = t0
where y and y' are the n vectors and the system consists of n equations. t0, y0 and y'0 are the given initial values of t, y and y', respectively.
This program is based on the Runge-Kutta-Nystroem method of order 4(3). See for details of coefficients (RKN4(3)4FM) and interpolant for dense output in the references (1) and (2), respectively.
This program is the reverse communication version of dopn43.
- Parameters
-
| [in] | n | Number of differential equations. (n >= 1) |
| [in,out] | t | Independent variable t. This program integrates from the initial value of t to tend. Depending on the setting of mode, this program may return at tout or every successful step to provide an intermediate result. [in] An initial value of t (t0).
[out] t = tend if the integration has been completed. t = tout or t of the latest step depending on mode in the case of intermediate return. |
| [in,out] | y[] | Array y[ly] (ly >= n)
Dependent variable y.
[in] Initial value of y (y0) at initial t (t0).
[out] The value of y (computed numerical solution) at t. |
| [in,out] | yp[] | Array yp[lyp] (lyp >= n)
[in] Initial values (y'0) of the derivatives y' at initial t (t0).
[out] The values of of the derivatives y' (computed numerical solution) at t. |
| [in] | tout | Set tout to the point (t) at which an intermediate result is desired in mode = 2 or 3. tout is not referenced in mode = 0 or 1.
In mode 2 or 3, y at tout is computed and returned with info = 2 or 3 respectively. To continue the integration to get result at new tout, it is possible to call this program again with new tout and without changing other variables including info.
t < tout <= tend is required. However, if the integration is in backward direction, tend <= tout < t is required. If tend is reached before tout, the integration is terminated immediately and t = tend and info = 0 is returned. |
| [in] | tend | The point at which an integration is completed. If tend is reached, the integration is terminated and t = tend and info = 0 is returned. The integration is possible for both forward (tend > t) and backward (tend < t) direction. |
| [in] | rtol | Scalar if itol = 0, or array rtol[lrtol] if itol = 1 (lrtol >= n) (rtol or rtol[i] >= 0)
The relative error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the itol.
rtol is used together with atol in a local error test at each step. The integration is performed using step sizes which are automatically selected so as to achieve the following equations.
If itol = 0 (i = 0 to n - 1): (local error) <= max(rtol*abs(y[i]) + atol, rtol*abs(yp[i]) + atol) (if iwork[11] = 0)
(local error) <= rtol*abs(y[i]) + atol (if iwork[11] = 1)
If itol = 1 (i = 0 to n - 1): (local error) <= max(rtol[i]*abs(y[i]) + atol, rtol[i]*abs(yp[i]) + atol) (if iwork[11] = 0)
(local error) <= rtol[i]*abs(y[i]) + atol[i] (if iwork[11] = 1)
Both rtol and atol, or, rtol[i] and atol[i] (i = 0 to n - 1) must not be 0 at the same time. |
| [in] | atol | Scalar if itol = 0, or array atol[latol] if itol = 1 (latol >= n) (atol or atol[i] >= 0)
The absolute error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the itol.
atol is used together with rtol in a local error test at each step (Refer to rtol above). |
| [in] | itol | Specifies whether the parameters rtol and atol are scalars or arrays.
= 0: rtol and atol are scalars (or arrays with one element).
= 1: rtol and atol are arrays. |
| [in] | mode | Mode of operation.
This program integrates from t0 to tend. Four modes of operation are provided depending on the timing of returning imtermediate results. When tend is reached, the integration is terminated and info = 0 is returned in any of four modes. = 0: Not return until tend. tout is not referenced.
= 1: Returns at every successful step (info = 1 is returned). tout is not referenced. The end point of the last step is returned in t. In the final step, the step size is adjusted to fit into tend.
= 2: Returns at tout during the ingegration to provide an intermediate result (t = tout and info = 2 are returned). To continue the integration to get result at new tout, it is possible to call this program again with new tout. In the final step, the step size is adjusted to fit into tout.
= 3: Returns at tout during the ingegration to provide an intermediate result (t = tout and info = 3 are returned). To continue the integration to get result at new tout, it is possible to call this program again with new tout. Different from mode = 2, the value of y at tout is computed by interpolation, and the step size is not adjusted even in the last step before tout. The steps through integration are same with those of mode = 1. |
| [in,out] | work[] | Array work[lwork]
Work array.
work[0] to work[19] serve as parameters for the program. If the input parameter is set to 0, the default parameter value defined for each parameter will be loaded.
[in]
work[0]: Initial step size (default = estimated by the code)
work[1]: Maximum step size (default = abs(tout - t))
work[4] and work[5]: Parameters for step size selection. The new step size is chosen subject to the restriction work[4] <= hnew/hold <= work[5]. (default work[4] = 0.2, work[5] = 10)
work[8]: The safety factor in step size prediction (0.0001 < work[8] < 1). (default = 0.9)
[out]
work[0]: Last step size. |
| [in] | lwork | Size of array work[]. (lwork >= 8*n + 40)
If lwork < 0, abs(lwork) will be used and work[0] to work[19] will be initialized to zeros. |
| [in,out] | iwork[] | Array iwork[liwork]
Work array.
iwork[0] to iwork[19] serve as parameters for the program. If the input parameter is set to 0, the default parameter value defined for each parameter will be loaded.
[in]
iwork[1]: Maximum number of allowed steps. (default = 10000)
iwork[11]: Specifies how determine the acceptance of the step. (default = 0)
= 0: Determine the acceptance of the step using the estimated errors of both y and y'.
= 1: Determine the acceptance of the step using the estimated error of y only.
iwork[12]: Specifies when counters iwork[13] to iwork[17] are reset to zero. (default = 0)
= 0: Reset if info = 0.
= 1: Do not reset even if info = 0.
[out]
Statistical information (counters).
iwork[13]: Number of function evaluations.
iwork[15]: Number of all computed steps.
iwork[16]: Number of accepted steps.
iwork[17]: Number of rejected steps. |
| [in] | liwork | Size of array iwork[]. (liwork >= 30)
If liwork < 0, abs(liwork) will be used and iwork[0] to iwork[19] will be initialized to zeros. |
| [in,out] | info | [in] Control code.
= 0: Set info = 0 on the initial call to start new problem. All variables will be initialized and the computation will begin.
= 1, 2, 3: When returned with info = 1, 2 or 3, it is possible to call this program again with new tout and without changing info to continue the integration.
[out] Return code.
= 0: Successful exit. Integration to tend completed.
< 0: The (-info)-th argument is invalid.
= 1: Returned to provide intermediate result in mode = 1. It is possible to call this program again to forward to the next step.
= 2: Returned at t = tout to provide intermediate result in mode = 2. It is possible to call this program again with new tout.
= 3: Returned at t = tout to provide intermediate result in mode = 3. It is possible to call this program again with new tout.
= 11: (Error) Maximum number of steps exceeded.
= 12: (Error) Step size becomes too small. |
| [out] | tt | irev = 1〜2: Refer to yypp[]. |
| [out] | yy[] | Array yy[lyy] (lyy >= n)
irev = 1〜2: Refer to yypp[]. |
| [in] | yypp[] | Array yypp[lyypp] (lyypp >= n)
irev = 1〜2: Computed second derivatives d2y/dt2 = f2(t, y) at given t = tt and y = yy[] should be set in yypp[] in the next call. |
| [in,out] | irev | Control variable for reverse communication interface.
[in] irev should be initialized to zero in the first call. In succeeding calls, irev should not be altered.
[out] If irev is not zero, complete the following tasks and call this program again without changing irev.
= 0: Computation finished. Check return code in info.
= 1〜2: Set the computed second derivatives at t = tt and y = yy[] in yypp[]. Do not alter any variables other than yypp[]. |
- Reference
- (1) J. R. Dormand, M. E. A. El-Mikkawy, P. J. Prince, "Families of Runge-Kutta-Nystrom Formulae", IMA J. of Numerical Analysis, 7, 235-250 (1987).
(2) J. R. Dormand, P. J. Prince "Runge-Kutta-Nystrom Triples", Comput. Math. Applic. Vol. 13, No. 12, 937-949 (1987).
|