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)