pracma (version 1.3.8)

cranknic: Crank-Nicolson Method

Description

The Crank-Nicolson method for solving ordinary differential equations is a combination of the generic steps of the forward and backward Euler methods.

Usage

cranknic(f, t0, t1, y0, ..., N = 100)

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, t1

start and end points of the interval.

y0

starting values as row or column vector; for $m$ equations y0 needs to be a vector of length m.

N

number of steps.

...

Additional parameters to be passed to the function.

Value

• 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.

Details

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.

References

Quarteroni, A., and F. Saleri (2006). Scientific Computing With MATLAB and Octave. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

ode23, newmark

Examples

##  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)

# 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()