diffEq (version 1.0-1)

stability.multistep: Plots the stability function of multistep methods

Description

...

Usage

stability.multistep (alpha, beta, add = FALSE, fill = NULL, ...) 
stability.bruteforce (Rez = seq(-2, 2, by = 0.02),  Imz = seq(-2, 2, by = 0.02),  func = function (z) return(abs(1 + h*z) <=1), fill =" "grey"," cex =" 1.5," add =" FALSE," ...)<="" div="">

Arguments

alpha

alpha coefficients of the multistep method.

beta

beta coefficients of the multistep method.

add

if TRUE, the new region will be added to the existing plot

fill

color of region to be filled

Rez

The range in the real plane for testing stability

Imz

The range in the imaginary plane for testing stability

func

The function to be tested; default is test for the euler method.

cex

The relative size of the plotting symbol. If too small, the region will not be completely covered. If too large, it will extend beyond its boundaries.

...

arguments passed to the plotting function.

Value

A matrix with the boundary value.

See Also

rkMethodPlot for plotting runge-kutta method steps.

Examples

par(mfrow=c(2,2))

## =============================================================================
## Stability regions for multistep methods
## =============================================================================

# Adams-Bashforth
stability.multistep(alpha = AdamsBashforth$alpha[2,], beta = AdamsBashforth$beta[2,], 
             xlim = c(-3,1), ylim = c(-1.5, 1.5),
             fill = "black",  main = "Adams-Bashforth")

stability.multistep(alpha = AdamsBashforth$alpha[3,], beta = AdamsBashforth$beta[3,], 
             add = TRUE, lty = 1, fill = "darkgrey")

stability.multistep(alpha = AdamsBashforth$alpha[4,], beta = AdamsBashforth$beta[4,],
             add = TRUE, fill = "lightgrey", lty = 1)

stability.multistep(alpha = AdamsBashforth$alpha[5,], beta = AdamsBashforth$beta[5,],
             add = TRUE, fill = "white", lty = 1)

legend("topleft", fill = c("black", "darkgrey", "lightgrey", "white"),
       title = "order", legend = 2:5)
writelabel("A")


# Adams-Moulton
stability.multistep(alpha = AdamsMoulton$alph[3,], beta = AdamsMoulton$beta[3,], 
             xlim = c(-8, 1), ylim = c(-4, 4), 
             fill = "black",   main = "Adams-Moulton")
stability.multistep(alpha = AdamsMoulton$alpha[4,], beta = AdamsMoulton$beta[4,], 
             add = TRUE, fill = "darkgrey")
stability.multistep(alpha = AdamsMoulton$alpha[5,], beta = AdamsMoulton$beta[5,], 
             add = TRUE, fill = "lightgrey")

legend("topleft", fill = c("black", "darkgrey", "lightgrey"),
       title = "order", legend = 3:5 )
writelabel("B")

# 2nd-order BDF
plot(0, type="n", xlim = c(-3, 12), ylim = c(-8, 8),
     main = "BDF order 2", 
     xlab = "Re(z)", ylab = "Im(z)")
rect(-100, -100, 100, 100, col = "lightgrey")
box()

stability.multistep (alpha = BDF$alpha[2,], beta = BDF$beta[2,], 
        fill = "white", add = TRUE)

writelabel("C")

# 4th-order BDF
plot(0, type="n", xlim=c(-3, 12), ylim = c(-8, 8),
     main = "BDF order 4", 
     xlab = "Re(z)", ylab = "Im(z)")
rect(-100, -100, 100, 100, col = "lightgrey")
box()
stability.multistep (alpha = BDF$alpha[4,], beta = BDF$beta[4,], 
        fill = "white", add = TRUE)

writelabel("D")

## =============================================================================
## Stability regions for runge-kutta methods
## =============================================================================
# Drawing the stability regions  - brute force
# stability function for explicit runge-kutta's
rkstabfunc <- function (z, order = 1) {
   h <- 1
   ss <- 1
   for (p in 1: order) ss <- ss + (h*z)^p / factorial(p)
   return (abs(ss) <= 1) 
 }

# regions for stability orders 5 to 1 - rather crude
Rez <- seq(-5, 1, by = 0.1)
Imz <- seq(-3, 3, by = 0.1) 

stability.bruteforce(main = "Explicit RK", 
   func = function(z) rkstabfunc(z, order = 5),
   Rez = Rez, Imz = Imz, fill = grey(0.95)  )

stability.bruteforce(add = TRUE, 
   func = function(z) rkstabfunc(z, order = 4),
   Rez = Rez, Imz = Imz, fill = grey(0.75)  )

stability.bruteforce(add = TRUE, 
   func = function(z) rkstabfunc(z, order = 3),
   Rez = Rez, Imz = Imz, fill = grey(0.55)  )

stability.bruteforce(add = TRUE, 
   func = function(z) rkstabfunc(z, order = 2),
   Rez = Rez, Imz = Imz, fill = grey(0.35)  )

stability.bruteforce(add = TRUE, 
   func = function(z) rkstabfunc(z, order = 1),
   Rez = Rez, Imz = Imz, fill = grey(0.15)  )

legend("topleft", legend = 5:1, title = "order", 
  fill = grey(c(0.95, 0.75, 0.55, 0.35, 0.15)))
  
# stability function for radau method

stability.bruteforce(main = "Radau 5", 
   Rez = seq(-5,15,by=0.1), Imz = seq(-10,10,by=0.12),
   func = function(z) return(abs((1 + 2*z/5 + z^2/20) /
                                 (1 - 3*z/5 + 3*z^2/20 - z^3/60)) <= 1),
   col = grey(0.8)  )