# R/odesolv.R In fda: Functional Data Analysis

#### Documented in odesolv

```odesolv <- function(bwtlist, ystart=diag(rep(1,norder)),
h0=width/100, hmin=width*1e-10, hmax=width*0.5,
EPS=1e-4, MAXSTP=1000)
{
#  Solve L u = 0,
#  L being an order M homogeneous linear differential operator,
#     (Lu)(t) = w_1(t) u(t) + w_2(t) Du(t) + ...
#                  w_m(t) D^{m-1} u(t) + D^m u(t) = 0
#  for function u and its derivatives up to order m - 1.
#  Each such solution is determined by the values of u and its
#  m - 1 derivatives at time t = 0.  These initial conditions are
#  contained in the columns of the matrix YSTART, which has exactly
#  m rows.  The number of solutions computed by ODESOLV is equal to
#  the number of columnsof YSTART.  In order for the solutions to be
#  linearly independent functions, the columns of YSTART must be
#  linearly independent.  This means that the maximum number of
#  linearly independent solutions is m.
#  The solution for each value of t is a matrix, y(t).
#  Any column of y(t) contains u, Du, ... , D^{m-1} u at argument
#  value t, for the corresponding set of starting values for these
#  first m derivatives. It is the job of this function to estimate
#  these values, and ODESOLV will choose a set of values of TNOW at
#  which these can be estimated with satisfactory accuracy.
#  ODESOLV uses the Runge-Kutta method, which is a good general
#  purpose numerical method.  But it does not work well for stiff
#  systems, and it can fail for poor choices of initial conditions
#  as well as other problems.

#  Arguments:
#  BWTLIST ... a list containing m = 1 weightfunctions.
#             The weight functions w_1, ... , w_m are the functions with
#             indices 2, 3, ..., m+1.
#  YSTART ... initial values for Y.  This is a matrix with M rows,
#             were M is the order of the operator L.  Any column of M
#             specifies intial values for derivatives 0, 1, ... M-1.  Each
#             column must specify a unique set of initial conditions.  A
#             frequent choice is the identity matrix of order M, and this is
#             the default.
#  H0     ... initial trial step size
#  HMIN   ... minimum step size
#  HMAX   ... maximum step size
#  EPS    ... error tolerance
#  MAXSTP ... maximum number of Runge Kutta steps permitted. If the equation
#             is difficult to solve, this may have to be increased.

#  Returns:
#  TP     ... vector of TRUE values used
#  YP     ... m by m by length(TP) array of Y-values generated for
#             values in TP.

#  Note that ODESOLV calls function DERIVS in order to evaluate the
#  differential operator.  Also, it works by redefining the order m
#  linear differential equation as a linear system of m first order
#  differential equations.  DERIVS evaluates the right side of this
#  linear system.

MAXWARN <- 10

#  determine the order of the system m

norder <- length(bwtlist)

#  determine the range of values over which the equation is solved

bfdPar1  <- bwtlist[[1]]
wbasis   <- bfdPar1\$fd\$basis
rangeval <- wbasis\$rangeval
tbeg     <- rangeval[1]
tend     <- rangeval[2]
width    <- tend - tbeg
tnow     <- tbeg

h        <- min(c(h0, hmax))
tp       <- tnow

#  set up the starting values

ystartd  <- dim(ystart)
if (ystartd[1] != norder) stop("YSTART has incorrect dimensions")
n        <- ystartd[2]
yp       <- c(ystart)
index    <- abs(yp) > 1e-10
y        <- ystart

#  Since ODESOLVE is slow, progress is displayed

cat("Progress:  each dot is 10% of interval\n")

#  initialize the solution

tdel  <- (tend - tbeg)/10
tdot  <- tdel
iwarn <- 0

#  solve the equation using a maximum of MAXSTP steps

for (nstp in 1:MAXSTP) {
#  evaluate the right side at the current value
dydt  <- derivs(tnow, y, bwtlist)
yscal <- c(abs(y) + abs(h*dydt) + 1e-30)[index]
if (nstp > 1) {
tp <- c(tp,tnow)
yp <- c(yp,c(y))
}
if (tnow >= tdot) {
cat(".")
tdot <- tdot + tdel
}
if ((tnow+h-tend)*(tnow+h-tbeg) > 0) h <- tend-tnow
#  take a Runge-Kutta step to the next value
result <-  rkqs(y, dydt, tnow, h, bwtlist, yscal, index, EPS)
tnow   <- result[[1]]
y      <- result[[2]]
h      <- result[[3]]
hnext  <- result[[4]]
#  test to see if interval has been covered, and return
#  if it has.
if ((tnow-tend)*(tend-tbeg) >= 0) {
cat(".")
tp <- c(tp,tnow)
yp <- c(yp,c(y))
yp <- array(yp,c(norder,n,length(tp)))
return(list(tp, yp))
}
#  test if the step is too small.
if (abs(hnext) < hmin) {
warning("Stepsize smaller than minimum")
hnext <- hmin
iwarn <- iwarn + 1
if (iwarn >= MAXWARN) stop("Too many warnings.")
}
#  test if the step is too large.
h <- min(c(hnext, hmax))
}
warning("Too many steps.")
}

#  ---------------------------------------------------------------

rkqs <- function(y, dydt, tnow, htry, bwtlist, yscal, index, EPS)
{
#  Take a single step using the Runge-Kutta procedure to
#  control for accuracy.  The function returns the new
#  value of t, the value of the solution at t, the step
#  size, and a proposal for the next step size.
h <- htry
#  check the accuracy of the step
result <- rkck(y, dydt, tnow, h, bwtlist)
ytemp  <- result[[1]]
yerr   <- c(result[[2]])[index]
errmax <- max(abs(yerr/yscal))/EPS
#  modify the step size if ERRMAX is too large
while (errmax > 1) {
htemp <- 0.9*h*(errmax^(-0.25))
h     <- max(c(abs(htemp),0.1*abs(h)))
tnew  <- tnow + h
if (tnew == tnow) stop("stepsize underflow in rkqs")
result <- rkck(y, dydt, tnow, h, bwtlist)
ytemp  <- result[[1]]
yerr   <- result[[2]]
errmax <- max(abs(yerr/yscal))/EPS
}
# set up the proposed next step
if (errmax > 1.89e-4) {
hnext <- 0.9*h*(errmax^(-0.2))
} else {
hnext <- 5.*h
}
tnow <- tnow + h
y <- ytemp
return( list(tnow, y, h, hnext) )
}

#  ------------------------------------------------------------------

rkck <- function(y, dydt, tnow, h, bwtlist)
{
#  Take a single Runge-Kutta step.
#  Return the solution, and an estimate of its error.
C1  <-   37/378
C3  <-  250/621
C4  <-  125/594
C6  <-  512/1771
DC5 <- -277/14336
DC1 <-  C1 - 2825/27648
DC3 <-  C3 - 18575/48384
DC4 <-  C4 - 13525/55296
DC6 <-  C6 - 0.25

ytemp <- y + h*0.2*dydt
ak2   <- derivs(tnow+0.2*h,   ytemp, bwtlist)
ytemp <- y + h*(0.075*dydt  + 0.225*ak2)
ak3   <- derivs(tnow+0.3*h,   ytemp, bwtlist)
ytemp <- y + h*(0.3*dydt    - 0.9*ak2 + 1.2*ak3)
ak4   <- derivs(tnow+0.6*h,   ytemp, bwtlist)
ytemp <- y + h*(-11*dydt/54 + 2.5*ak2 + (-70*ak3+35*ak4)/27)
ak5   <- derivs(tnow+h,       ytemp, bwtlist)
ytemp <- y + h*(1631*dydt/55296  + 575*ak2/512 + 575*ak3/13824 +
44275*ak4/110592 + 253*ak5/4096)
ak6   <- derivs(tnow+0.875*h, ytemp, bwtlist)
yout  <- y + h*(C1*dydt + C3*ak3 + C4*ak4 + C6*ak6)
yerr  <- h*(DC1*dydt    + DC3*ak3 + DC4*ak4 + DC5*ak5 + DC6*ak6)
return (list( yout, yerr ) )
}
```

## Try the fda package in your browser

Any scripts or data that you put into this service are public.

fda documentation built on May 29, 2024, 11:26 a.m.