ode23 | R Documentation |
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.
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)
f |
function in the differential equation |
t0 , tfinal |
start and end points of the interval. |
y0 |
starting values as column vector;
for |
jac |
jacobian of |
rtol , atol |
relative and absolute tolerance. |
hmax |
maximal step size, default is |
... |
Additional parameters to be passed to the function. |
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.
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.
Copyright (c) 2004 C. Moler for the Matlab textbook version ode23tx
.
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. https://www.mathworks.com/moler/chapters.html.
rk4sys
, deval
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.