Description Usage Arguments Value Author(s) See Also Examples
...
1 2 3 4 5 6 |
alpha |
alpha coefficients of the multistep method. |
beta |
beta coefficients of the multistep method. |
add |
if |
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 |
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. |
A matrix with the boundary value.
Karline Soetaert
rkMethodPlot for plotting runge-kutta method steps.
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) )
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.