trapezoidal <- function ( f, t0, tstop, y0, n, ... ) #*****************************************************************************80 # ## trapezoidal() implements the trapezoidal ODE solver. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 April 2021 # # Author: # # Adapted from "cranknic.R, Hans Borchers, the pracma library. # Modifications by John Burkardt # # Input: # # function f ( t, y, ... ): evaluates the right hand side of the ODE. # # real T0: the initial time. # # real TSTOP: the final time. # # real Y0: the initial solution values. # # integer N: the number of steps to take. # # ...: optional additional arguments to be passed to F(). # # Output: # # list ( T(1:n+1), Y(1:n+1) ): the computed sequence of # solution estimates. # { stopifnot ( is.numeric(y0), is.numeric(t0), length(t0) == 1, is.numeric(tstop), length(tstop) == 1 ) if ( is.vector(y0) ) { y0 <- as.matrix(y0) } else if ( is.matrix(y0) ) { if ( ncol(y0) != 1 ) { stop ( "trapezoidal: Argument 'y0' must be a row or column vector.") } } fun <- match.fun(f) f <- function(t, y) fun ( t, y, ... ) y_length <- length ( y0 ) y <- y0 yout <- matrix ( NA, n, y_length ) yout[1, ] <- c(y0) dt <- ( tstop - t0 ) / ( n - 1 ) t <- 0 ts <- linspace ( t0, tstop, n ) # # internal function used for root finding # cnfun <- function(w) w - y - 0.5* dt * ( f(t, w) + f(t, y) ) f_length <- length ( f ( t0, y0 ) ) if ( f_length != y_length ) { stop ( "trapezoidal: Function f must return a vector the same length as 'y0'.") } # # Choose the solver used for root finding. # if ( y_length == 1 ) { solver <- fzero } else { solver <- fsolve } for ( i in 2 : n ) { t <- ts[i] w <- solver ( cnfun, y )\$x yout[i, ] <- w y <- w } if ( y_length == 1 ) { yout <- drop ( yout ) } return ( list ( t = ts, y = yout ) ) }