Stability: Plots the stability function of multistep methods

Description Usage Arguments Value Author(s) See Also Examples

Description

...

Usage

1
2
3
4
5
6
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, ...)   

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.

Author(s)

Karline Soetaert

See Also

rkMethodPlot for plotting runge-kutta method steps.

Examples

  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
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
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)  )

diffEq documentation built on May 2, 2019, 2:18 a.m.