PDE imaginary
## =============================================================================
## 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 partial differential equation model with complex numbers
library(deSolve)
alf <- 0.5
gam <- 1
Schrodinger <- function(t, u, parms) {
du <- 1i * tran.1D (C = u, D = 1, dx = xgrid)$dC +
1i * gam * abs(u)^2 * u
list(du)
}
N <- 300
xgrid <- setup.grid.1D(-20, 80, N = N)
x <- xgrid$x.mid
c1 <- 1
c2 <- 0.1
sech <- function(x) 2/(exp(x) + exp(-x))
soliton <- function (x, c1)
sqrt(2*alf/gam) * exp(0.5*1i*c1*x) * sech(sqrt(alf)*x)
yini <- soliton(x, c1) + soliton(x-25, c2)
times <- seq(0, 40, by = 0.1)
print(system.time(
out <- ode.1D(y = yini, parms = NULL, func = Schrodinger,
times = times, dimens = 300, method = "adams")
))
image(abs(out), grid = x, ylab = "x", main = "two solitons",
legend = TRUE)
Go back