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

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.

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

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.

Examples

Run this code

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

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning