cranknic | R Documentation |
The Crank-Nicolson method for solving ordinary differential equations is a combination of the generic steps of the forward and backward Euler methods.
cranknic(f, t0, t1, y0, ..., N = 100)
f |
function in the differential equation |
t0 , t1 |
start and end points of the interval. |
y0 |
starting values as row or column vector;
for |
N |
number of steps. |
... |
Additional parameters to be passed to the function. |
Adding together forward and backword Euler method in the cranknic
method is by finding the root of the function merging these two formulas.
No attempt is made to catch any errors in the root finding functions.
List with components t
for grid (or ‘time’) points between t0
and t1
, and y
an n-by-m matrix with solution variables in
columns, i.e. each row contains one time stamp.
This is for demonstration purposes only; for real problems or applications
please use ode23
or rkf54
.
Quarteroni, A., and F. Saleri (2006). Scientific Computing With MATLAB and Octave. Second Edition, Springer-Verlag, Berlin Heidelberg.
ode23
, newmark
## Newton's example
f <- function(x, y) 1 - 3*x + y + x^2 + x*y
sol100 <- cranknic(f, 0, 1, 0, N = 100)
sol1000 <- cranknic(f, 0, 1, 0, N = 1000)
## Not run:
# Euler's forward approach
feuler <- function(f, t0, t1, y0, n) {
h <- (t1 - t0)/n; x <- seq(t0, t1, by = h)
y <- numeric(n+1); y[1] <- y0
for (i in 1:n) y[i+1] <- y[i] + h * f(x[i], y[i])
return(list(x = x, y = y))
}
solode <- ode23(f, 0, 1, 0)
soleul <- feuler(f, 0, 1, 0, 100)
plot(soleul$x, soleul$y, type = "l", col = "blue",
xlab = "", ylab = "", main = "Newton's example")
lines(solode$t, solode$y, col = "gray", lwd = 3)
lines(sol100$t, sol100$y, col = "red")
lines(sol1000$t, sol1000$y, col = "green")
grid()
## System of differential equations
# "Herr und Hund"
fhh <- function(x, y) {
y1 <- y[1]; y2 <- y[2]
s <- sqrt(y1^2 + y2^2)
dy1 <- 0.5 - 0.5*y1/s
dy2 <- -0.5*y2/s
return(c(dy1, dy2))
}
sol <- cranknic(fhh, 0, 60, c(0, 10))
plot(sol$y[, 1], sol$y[, 2], type = "l", col = "blue",
xlab = "", ylab = "", main = '"Herr und Hund"')
grid()
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.