PDE
## =============================================================================
## R-Code to solve example 9.3.4 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.
## Solving differential equations in R. UseR, Springer, 248 pp.
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
## implemented by Karline Soetaert
## =============================================================================
## A 3-D partial differential equation model
library(deSolve)
Nx <- 100
Ny <- 100
xgrid <- setup.grid.1D(-7, 7, N=Nx)
ygrid <- setup.grid.1D(-7, 7, N=Ny)
x <- xgrid$x.mid
y <- ygrid$x.mid
sinegordon2D <- function(t, C, parms) {
u <- matrix(nrow = Nx, ncol = Ny,
data = C[1 : (Nx*Ny)])
v <- matrix(nrow = Nx, ncol = Ny,
data = C[(Nx*Ny+1) : (2*Nx*Ny)])
dv <- tran.2D (C = u, C.x.up = 0, C.x.down = 0,
C.y.up = 0, C.y.down = 0,
D.x = 1, D.y = 1,
dx = xgrid, dy = ygrid)$dC - sin(u)
list(c(v, dv))
}
peak <- function (x, y, x0 = 0, y0 = 0)
exp(-((x-x0)^2 + (y-y0)^2))
uini <- outer(x, y,
FUN = function(x, y) peak(x, y, 2,2) + peak(x, y,-2,-2)
+ peak(x, y,-2,2) + peak(x, y, 2,-2))
vini <- rep(0, Nx*Ny)
times <- 0:3
print(system.time(
out <- ode.2D (y = c(uini, vini), times = times,
parms = NULL, func = sinegordon2D,
names = c("u", "v"),
dimens = c(Nx, Ny), method = "ode45")
))
mr <- par(mar = c(0, 0, 1, 0))
image(out, main = paste("time =", times), which = "u",
grid = list(x = x, y = y), method = "persp",
border = NA, col = "grey", box = FALSE,
shade = 0.5, theta = 30, phi = 60, mfrow = c(2, 2),
ask = FALSE)
par(mar = mr)
Go back