Description Usage Author(s) See Also Examples
Coefficients alpha
and beta
of the Adams-Bashforth,
Adams-Moulton and Backward differentiation formulae.
sum_j=0^k alpha_j y_n-j = h sum_(j=0)^k beta_j f(x_n-j,y_n-j)
For the BDF methods, the angle
of the stability-region (the alpha
of the A(alpha) stability, in radians is also given.
1 2 3 |
Karline Soetaert
stability.multistep for plotting stability regions
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 | ## =============================================================================
## Stability properties
## =============================================================================
BDF
stability.multistep(alpha = BDF$alpha[3,], beta = BDF$beta[3,],
xlim = c(-7,7), ylim = c(-7,7))
stability.multistep(alpha = AdamsMoulton$alpha[3,], beta = AdamsMoulton$beta[3,],
xlim = c(-7,7), ylim = c(-7,7) )
stability.multistep(alpha = AdamsBashforth$alpha[3,], beta = AdamsBashforth$beta[3,] )
## =============================================================================
## Running a BDF
## =============================================================================
# test model
ode1 <- function (t, y) return(cos(t)*y )
h <- 0.01
times <- seq(from = 0, to = 20, by = h)
yout <- vector (length = length(times))
yout[1] <- 1
# 3rd order BDF
Alpha <- BDF$alpha [3,2:4]
Beta <- BDF$beta[3,]
bdf <- function(y, t, h, f, ys) {
rootfun <- function(ynext)
- ynext - sum(Alpha * ys) + Beta * h * f(t + h, ynext)
y <- multiroot(f=rootfun, start=y)$root
ys[2:3] <- ys[1:2]
ys[1] <- y
list (y = y, ys=ys)
}
# Spinup uses Euler...
Euler <- function(y, t, h, f) {
fn <- f(t, y)
ynext <- y + h * fn
list (y = ynext, fn = fn)
}
for (i in 2:3)
yout[i] <- Euler(yout[i-1], times[i-1], h, ode1)$y
ys <- rev(yout[1:3])
# BDF steps
for (i in 4:length(times)){
step <- bdf (y=yout[i-1], t=times[i-1], h=h, f=ode1, ys=ys)
yout[i] <- step$y
ys <- step$ys
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.