Learn R Programming

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

Run this code
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)  )

Run the code above in your browser using DataLab