# Stability: Plots the stability function of multistep methods In diffEq: Functions from the book Solving Differential Equations in R

...

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

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.