# Non-stiff (and stiff) ODE solvers Runge-Kutta (2, 3)-method with variable step size, resp. (4,5)-method with Dormand-Price coefficients, or (7,8)-pairs with Fehlberg coefficients. The function `f(t, y)` has to return the derivative as a column vector. ```r ode23(f, t0, tfinal, y0, ..., rtol = 1e-3, atol = 1e-6) ode23s(f, t0, tfinal, y0, jac = NULL, ..., rtol = 1e-03, atol = 1e-06, hmax = 0.0) ode45(f, t0, tfinal, y0, ..., atol = 1e-6, hmax = 0.0) ode78(f, t0, tfinal, y0, ..., atol = 1e-6, hmax = 0.0) ``` ## Arguments - `f`: function in the differential equation $y' = f(x, y)$; defined as a function $R \times R^m \rightarrow R^m$, where $m$ is the number of equations. - `t0, tfinal`: start and end points of the interval. - `y0`: starting values as column vector; for $m$ equations `u0` needs to be a vector of length `m`. - `jac`: jacobian of `f` as a function of `x` alone; if not specified, a finite difference approximation will be used. - `rtol, atol`: relative and absolute tolerance. - `hmax`: maximal step size, default is `(tfinal - t0)/10.` - `...`: Additional parameters to be passed to the function. ## Details `ode23` is an integration method for systems of ordinary differential equations using second and third order Runge-Kutta-Fehlberg formulas with automatic step-size. `ode23s` can be used to solve a stiff system of ordinary differential equations, based on a modified Rosenbrock triple method of order (2,3); See section 4.1 in [Shampine and Reichelt]. `ode45` implements Dormand-Prince (4,5) pair that minimizes the local truncation error in the 5th-order estimate which is what is used to step forward (local extrapolation). Generally it produces more accurate results and costs roughly the same computationally. `ode78` implements Fehlberg's (7,8) pair and is a 7th-order accurate integrator therefore the local error normally expected is O(h^8). However, because this particular implementation uses the 8th-order estimate for xout (i.e. local extrapolation) moving forward with the 8th-order estimate will yield errors on the order of O(h^9). It requires 13 function evaluations per integration step. ## Returns List with components `t` for grid (or `time') points between `t0` and `tfinal`, and `y` an n-by-m matrix with solution variables in columns, i.e. each row contains one time stamp. ## References Ascher, U. M., and L. R. Petzold (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. SIAM. L.F. Shampine and M.W. Reichelt (1997). The MATLAB ODE Suite. SIAM Journal on Scientific Computing, Vol. 18, pp. 1-22. Moler, C. (2004). Numerical Computing with Matlab. Revised Reprint, SIAM. See www.mathworks.com/moler/chapters.html. ## Note Copyright (c) 2004 C. Moler for the Matlab textbook version `ode23tx`. ## See Also `rk4sys`, `deval` ## Examples ```r ## Example1: Three-body problem f <- function(t, y) as.matrix(c(y[2]*y[3], -y[1]*y[3], 0.51*y[1]*y[2])) y0 <- as.matrix(c(0, 1, 1)) t0 <- 0; tf <- 20 sol <- ode23(f, t0, tf, y0, rtol=1e-5, atol=1e-10) ## Not run: matplot(sol$t, sol$y, type = "l", lty = 1, lwd = c(2, 1, 1), col = c("darkred", "darkblue", "darkgreen"), xlab = "Time [min]", ylab= "", main = "Three-body Problem") grid() ## End(Not run) ## Example2: Van der Pol Equation # x'' + (x^2 - 1) x' + x = 0 f <- function(t, x) as.matrix(c(x[1] * (1 - x[2]^2) -x[2], x[1])) t0 <- 0; tf <- 20 x0 <- as.matrix(c(0, 0.25)) sol <- ode23(f, t0, tf, x0) ## Not run: plot(c(0, 20), c(-3, 3), type = "n", xlab = "Time", ylab = "", main = "Van der Pol Equation") lines(sol$t, sol$y[, 1], col = "blue") lines(sol$t, sol$y[, 2], col = "darkgreen") grid() ## End(Not run) ## Example3: Van der Pol as stiff equation vdP <- function(t,y) as.matrix(c(y[2], 10*(1-y[1]^2)*y[2]-y[1])) ajax <- function(t, y) matrix(c(0, 1, -20*y[1]*y[2]-1, 10*(1-y[1]^2)), 2,2, byrow = TRUE) sol <- ode23s(vdP, t0, tf, c(2, 0), jac = ajax, hmax = 1.0) ## Not run: plot(sol$t, sol$y[, 1], col = "blue") lines(sol$t, sol$y[, 1], col = "blue") lines(sol$t, sol$y[, 2]/8, col = "red", lwd = 2) grid() ## End(Not run) ## Example4: pendulum m = 1.0; l = 1.0 # [kg] resp. [m] g = 9.81; b = 0.7 # [m/s^2] resp. [N s/m] fp = function(t, x) c( x[2] , 1/(1/3*m*l^2)*(-b*x[2]-m*g*l/2*sin(x[1])) ) t0 <- 0.0; tf <- 5.0; hmax = 0.1 y0 = c(30*pi/180, 0.0) sol = ode45(fp, t0, tf, y0, hmax = 0.1) ## Not run: matplot(sol$t, sol$y, type = "l", lty = 1) grid() ## End(Not run) ## Example: enforced pendulum g <- 9.81 L <- 1.0; Y <- 0.25; w <- 2.5 f <- function(t, y) { as.matrix(c(y[2], -g/L * sin(y[1]) + w^2/L * Y * cos(y[1]) * sin(w*t))) } y0 <- as.matrix(c(0, 0)) sol <- ode78(f, 0.0, 60.0, y0, hmax = 0.05) ## Not run: plot(sol$t, sol$y[, 1], type="l", col="blue") grid() ## End(Not run) ```

ode function