arf_rk: Numerical Solution of Systems of Ordinary Differential Equations

Description

Solves numerically the initial value problem $$y'(t) = f(t, y(t)),\quad y(0) = y_{0},$$ using an explicit, adaptive or non-adaptive Runge-Kutta method, by default the Dormand-Prince method.

The main difference between this function and function rk in package deSolve is the use of arf instead of double, enabling computation of times, solution values, and solution derivatives in arbitrary precision. A corollary is that users can choose arbitrarily small rtol and atol provided that the working precision is sufficiently high.

Interpolation is not yet implemented.

Usage

arf_rk(func, t, y0, param = NULL, rtol = NULL, atol = NULL,
       hmin = 0, hmax = Inf, hini = NULL, smax = NULL,
       method = .rk.method.dormand.prince(), progress = 0L,
       prec = flintPrec())

Value

A list with components:

t: an increasing arf vector storing time points.
y: an arf matrix with length(t) rows and length(y0) columns storing the numerical solution. In the event of early termination, trailing rows are filled with NaN.

Arguments

func: a function of the form function (t, y, param, prec) specifying the system. Unused trailing arguments can omitted.
t: an increasing numeric or arf vector storing time points.
y0: a numeric or arf vector storing the initial value.
param: an R object typically specifying parameters of the system, passed to func.
rtol: a positive number less than 1 controlling external step adaptation in adaptive methods. 2^(-prec/2) by default, unused by non-adaptive methods.
atol: a non-negative number less than 1 controlling external step adaptation in in adaptive methods. 2^(-prec/2) by default, unused by non-adaptive methods.
hmin: a non-negative number indicating a minimum external step size in adaptive methods. Early termination results if it seems that a smaller step size is needed to achieve sufficiently small error. The default value is 0, indicating that the step size can become arbitrarily small. Unused by non-adaptive methods.
hmax: a positive number indicating a maximum external step size in adaptive methods. The default value is Inf, indicating that the step size is bounded above by diff(t) and nothing else. Unused by non-adaptive methods.
hini: a positive number indicating the initial external step size in adaptive methods and the fixed external step size in non-adaptive methods. min(diff(t)) by default.
smax: a non-negative integer indicating a maximum number of internal steps per external step. Early termination results after smax * (length(t) - 1) internal steps. 256 * prec by default.
method: a list specifying a Runge-Kutta method.
progress: an integer flag determining how progress is indicated. . is printed after each external step if progress >= 1; o and x are printed after each accepted and rejected internal step if progress >= 2.
prec: a positive integer indicating the working precision as a number of bits, passed to func.

Details

func(t, y, param, prec) computes the derivative of the solution at time t given the value y of the solution at time t and optional parameters param, where t is an arf vector of length 1 and y is (and the return value of func must be) an arf vector of length equal to length(y0).

The list method must have exactly the following components for a d-stage method:

a: a numeric or fmpq or arf vector of length d*(d-1)/2 storing coefficients.
b, bb: numeric or fmpq or arf vectors of length d storing lower and higher order weights, each summing to 1. Set bb to NULL to specify a non-adaptive method.
c: a numeric or fmpq or arf vector of length d storing nodes for internal steps. The first element must be 0.
p: a positive integer giving the order of the method, such that the global error is $O(h^{p})$.

Examples

Run this code

F.linexp <- function (t, y) c(arf(1), -y[2L])
tt <- 0:10
y0 <- c(u = 1, v = 1)
L <- arf_rk(F.linexp, tt, y0)
L. <- list(t = tt, y = cbind(u = y0[1] + tt, v = y0[2] * exp(-tt)))
stopifnot(all.equal(L, L., check.class = FALSE))