https://github.com/cran/pracma
Raw File
Tip revision: 26e049d70b4a1c237987e260cba68f6a9413736c authored by Hans W. Borchers on 09 April 2019, 04:10:07 UTC
version 2.2.5
Tip revision: 26e049d
euler_heun.R
##
##  e u l e r H e u n . R  Euler-Heun ODE Solver
##


euler_heun <- function(f, a, b, y0, n = 100, improved = TRUE, ...) {
    stopifnot(is.numeric(a), is.numeric(b), length(a) == 1, length(b) == 1,
              is.numeric(y0), length(y0) == 1)

    fun <- match.fun(f)
    f <- function(t, y) fun(t, y, ...)
    if (length(f(a, y0)) != 1)
        stop("Argument function 'f' must be an univariate function.")

    h <- (b - a)/n
    t <- seq(a, b, length.out = n+1)
    y <- numeric(n+1)
    y[1] <- y0

    for (i in 1:n) {
        y[i+1] <- y[i] + h*f(t[i], y[i])
        if (improved) {
            y[i+1] <- y[i] + h * (f(t[i], y[i]) + f(t[i+1], y[i+1]))/2.0
        }
    }
    return(list(t = t, y = y))
}
back to top