From 400f851c026b1e18510463f393712607fdc21850 Mon Sep 17 00:00:00 2001 From: Burak Bayramli Date: Thu, 29 Aug 2019 16:13:00 +0300 Subject: [PATCH] - --- .../BVPinR.R | 313 ++++++++++++++++ .../DAEinR.R | 212 +++++++++++ .../DDEinR.R | 179 +++++++++ .../ODEinR.R | 340 +++++++++++++++++ .../PDEinR.R | 352 ++++++++++++++++++ .../light.rda | Bin 0 -> 1719 bytes 6 files changed, 1396 insertions(+) create mode 100644 Solving_Differential_Equations_in_R_Cash/BVPinR.R create mode 100644 Solving_Differential_Equations_in_R_Cash/DAEinR.R create mode 100644 Solving_Differential_Equations_in_R_Cash/DDEinR.R create mode 100644 Solving_Differential_Equations_in_R_Cash/ODEinR.R create mode 100644 Solving_Differential_Equations_in_R_Cash/PDEinR.R create mode 100644 Solving_Differential_Equations_in_R_Cash/light.rda diff --git a/Solving_Differential_Equations_in_R_Cash/BVPinR.R b/Solving_Differential_Equations_in_R_Cash/BVPinR.R new file mode 100644 index 000000000..b240fa09a --- /dev/null +++ b/Solving_Differential_Equations_in_R_Cash/BVPinR.R @@ -0,0 +1,313 @@ +## ============================================================================= +## +## Examples from the book +## Soetaert, K., Cash, J.R. and Mazzia, F. (2012). +## Solving Differential Equations in R. +## Springer +## +## Chapter 12. Solving Boundary Value problems in R. +## +## ============================================================================= + +## ----------------------------------------------------------------------------- +## A simple example +## ----------------------------------------------------------------------------- +library(bvpSolve) +prob7 <- function(x, y, pars) { + list(c( y[2], + 1/eps * (-x*y[2] + y[1] - (1+eps*pi*pi)* + cos(pi*x) - pi*x*sin(pi*x)))) +} + + +eps <- 0.1 +sol <- bvptwp(yini = c(y = -1, y1 = NA), + yend = c(1, NA), func = prob7, + x = seq(-1, 1, by = 0.01)) + +prob7_2 <- function(x, y, pars) { + list(1/eps * (-x*y[2] + y[1] - (1+eps*pi*pi)* + cos(pi*x) - pi*x*sin(pi*x))) +} +sol1 <- bvptwp(yini = c(y = -1, y1 = NA), + yend = c(1, NA), func = prob7_2, + order = 2, x = seq(-1, 1, by = 0.01)) +head(sol, n=3) +eps <-0.0005 +sol2 <- bvptwp(yini = c(y = -1, y1 = NA), + yend = c(1, NA), func = prob7, + x = seq(-1, 1, by=0.01)) + +plot(sol, sol2, col = "black", lty = c("solid", "dashed"), + lwd = 2) + +## ----------------------------------------------------------------------------- +## The swirling flow +## ----------------------------------------------------------------------------- +swirl <- function (t, Y, eps) { + with(as.list(Y), + list(c((g*f1 - f*g1)/eps, + (-f*f3 - g*g1)/eps)) + ) +} + +eps <- 0.001 +x <- seq(from = 0, to = 1, length = 100) +yini <- c(g = -1, g1 = NA, f = 0, f1 = 0, f2 = NA, f3 = NA) +yend <- c(1, NA, 0, 0, NA, NA) +Soltwp <- bvptwp(x = x, func = swirl, order = c(2, 4), + par = eps, yini = yini, yend = yend) + +pairs(Soltwp, main = "swirling flow III, eps=0.001") + +diagnostics(Soltwp) + +## ----------------------------------------------------------------------------- +## The musn problem +## ----------------------------------------------------------------------------- +musn <- function(x, Y, pars) { + with (as.list(Y), { + du <- 0.5 * u * (w - u) / v + dv <- -0.5 * (w - u) + dw <- (0.9 - 1000 * (w - y) - 0.5 * w * (w - u))/z + dz <- 0.5 * (w - u) + dy <- -100 * (y - w) + return(list(c(du, dv, dw, dz, dy))) + }) +} + +bound <- function(i, Y, pars) { + with (as.list(Y), { + if (i == 1) return (u - 1) + if (i == 2) return (v - 1) + if (i == 3) return (w - 1) + if (i == 4) return (z + 10) + if (i == 5) return (w - y) + }) +} + +xguess <- seq(0, 1, length.out = 5) +yguess <- matrix(ncol = 5, + data = (rep(c(1, 1, 1, -10, 0.91), 5))) +rownames(yguess) <- c("u", "v", "w", "z", "y") +xguess +yguess + +Sol <- bvptwp(x = x, func = musn, bound = bound, + xguess = xguess, yguess = yguess, + leftbc = 4, atol = 1e-10) + +yguess <- matrix(ncol = 5, data = (rep(c(1,1,1, 10, 0.91), 5))) +rownames(yguess) <- c("u", "v", "w", "z", "y") +Sol2<- bvpcol(x = x, func = musn, bound = bound, + xguess = xguess, yguess = yguess, + leftbc = 4, atol = 1e-10) + +plot(Sol, Sol2, which = "y", lwd = 2) + +## ----------------------------------------------------------------------------- +## The problem 19 +## ----------------------------------------------------------------------------- +Prob19 <- function(x, y, eps) { + pix = pi*x + list(c(y[2], + (pi/2*sin(pix/2)*exp(2*y[1])-exp(y[1])*y[2])/eps)) +} +x <- seq(0, 1, by = 0.01) +eps <- 1e-2 +mod1 <- bvptwp(func = Prob19, yini = c(0, NA), yend = c(0, NA), + x = x, par = eps) +diagnostics(mod1) + +plot(mod1, lwd = 2) + +## ----------------------------------------------------------------------------- +## The problem 19 (part b) +## ----------------------------------------------------------------------------- +xguess <- mod1[,1] +yguess <- t(mod1[,2:3]) + +eps <- 1e-3 +mod2 <- bvpcol(func = Prob19, yini = c(0, NA), yend = c(0, NA), + x = x, par = eps, + xguess = xguess, yguess = yguess) + +eps <- 1e-7 +mod2 <- bvpcol(func = Prob19, yini = c(0, NA), yend = c(0, NA), + x = x, par = eps, eps = eps, atol = 1e-4) + +diagnostics(mod2) + +## ----------------------------------------------------------------------------- +## The elastica problem +## ----------------------------------------------------------------------------- + +Elastica <- function (x, y, pars) { + list( c(cos(y[3]), + sin(y[3]), + y[4], + y[5] * cos(y[3]), + 0)) +} +bvpsol <- bvpcol(func = Elastica, + yini = c(x = 0, y = 0, p = NA, k = 0, F = NA), + yend = c(x = NA, y = 0, p = -pi/2, k = NA, F = NA), + x = seq(from = 0, to = 0.5, by = 0.01)) +bvpsol[1,"F"] + +plot(bvpsol, lwd = 2) + +## ----------------------------------------------------------------------------- +## The measel problem +## ----------------------------------------------------------------------------- + +measel <- function(t, y, pars) { + bet <- 1575 * (1 + cos(2 * pi * t)) + dy1 <- mu - bet * y[1] * y[3] + dy2 <- bet * y[1] * y[3] - y[2] / lam + dy3 <- y[2] / lam - y[3] / eta + dy4 <- 0 + dy5 <- 0 + dy6 <- 0 + list(c(dy1, dy2, dy3, dy4, dy5, dy6)) +} + +bound <- function(i, y, pars) { + if ( i == 1 | i == 4) return( y[1] - y[4]) + if ( i == 2 | i == 5) return( y[2] - y[5]) + if ( i == 3 | i == 6) return( y[3] - y[6]) +} + +mu <- 0.02 ; lam <- 0.0279 ; eta <- 0.1 +x <- seq(from = 0, to = 1, by = 0.01) + +Sola <- bvpshoot(func = measel, bound = bound, + x = x, leftbc = 3, atol = 1e-12, rtol = 1e-12, + guess = c(y1 = 1, y2 = 1, y3 = 1, y4 = 1, y5 = 1, y6 = 1)) + +yguess <- matrix(ncol = length(x), nrow = 6, data = 1) +rownames(yguess) <- paste("y", 1:6, sep="") +Sol <- bvptwp (func = measel, bound = bound, + x = x, leftbc = 3, xguess = x, yguess = yguess) +max(abs(Sol[1,-1] - Sol[nrow(Sol),-1])) + +plot(Sol, lwd = 2, which = 1:3, mfrow = c(1, 3)) + +## ----------------------------------------------------------------------------- +## The nerve impulse problem +## ----------------------------------------------------------------------------- + +nerve <- function (x, y, p) + list(c(3 * y[3] * (y[1] + y[2] - 1/3 * (y[1]^3) -1.3), + (-1/3) * y[3] * (y[1] - 0.7 + 0.8 * y[2]) , + 0, + 0, + 0) + ) + +bound <- function(i, y, p) { + if (i ==1) return (-y[3]* (y[1] - 0.7 + 0.8*y[2])/3 -1) + if (i ==2) return (y[1] - y[4] ) + if (i ==3) return (y[2] - y[5] ) + if (i ==4) return (y[1] - y[4] ) + if (i ==5) return (y[2] - y[5] ) +} +xguess <- seq(0, 1, by = 0.1) +yguess <- matrix(nrow = 5, ncol = length(xguess), data = 5.) +yguess[1,] <- sin(2 * pi * xguess) +yguess[2,] <- cos(2 * pi * xguess) +rownames(yguess) <- c("y1", "y2", "T", "y1ini", "y2ini") + +Sol <- bvptwp(func = nerve, bound = bound, + x = seq(0, 1, by = 0.01),leftbc = 3, + xguess = xguess, yguess = yguess) +Sol[1,] + +plot(Sol, lwd = 2, which = c("y1", "y2")) + +## ----------------------------------------------------------------------------- +## Integral constraints +## ----------------------------------------------------------------------------- +library(bvpSolve) + +integro <- function (t, u, p) + list(c(u[2], -p*exp(u[1]), u[2])) + +yini <- c(u1 = 1, u2 = NA, I = 0) +yend <- c(NA, NA, 1) +x <- seq(from = 0, to = 1, by = 0.01) + +out <- bvpcol (yini = yini, yend = yend, func = integro, + x = x, parms = 0.5) + +## ----------------------------------------------------------------------------- +## Sturm-Liouville +## ----------------------------------------------------------------------------- + +Sturm <- function(x, y, p) { + dy1 <- y[2] + dy2 <- -y[3] * y[1] + dy3 <- 0. + list( c(dy1, dy2, dy3)) +} + +yini <- c(y = 0, dy = 1, lambda = NA) +yend <- c(y = 0, dy = NA, lambda = NA) +x <- seq(from = 0, to = pi, by = pi/10) + +S1 <- bvpshoot(yini = yini, yend = yend, func = Sturm, + parms = 0, x = x) + +(lambda1 <- S1[1, "lambda"]) + +ana <- function(x, lambda) sin(x*sqrt(lambda))/sqrt(lambda) +max (abs(S1[,2]-ana(S1[,1],lambda1))) + +## ----------------------------------------------------------------------------- +## Reaction - transport model +## ----------------------------------------------------------------------------- +library(ReacTran) +N <- 1000 +Grid <- setup.grid.1D(N = N, L = 100000) +v <- 1000; D <- 1e7; O2s <- 300; +NH3in <- 500; O2in <- 100; NO3in <- 50 +r <- 0.1; k <- 1.; p <- 0.1 + +Estuary <- function(t, y, parms) { + NH3 <- y[1:N] + NO3 <- y[(N+1):(2*N)] + O2 <- y[(2*N+1):(3*N)] + tranNH3<- tran.1D (C = NH3, D = D, v = v, + C.up = NH3in, C.down = 10, dx = Grid)$dC + tranNO3<- tran.1D (C = NO3, D = D, v = v, + C.up = NO3in, C.down = 30, dx = Grid)$dC + tranO2 <- tran.1D (C = O2 , D = D, v = v, + C.up = O2in, C.down = 250, dx = Grid)$dC + + reaeration <- p * (O2s - O2) + r_nit <- r * O2 / (O2 + k) * NH3 + + dNH3 <- tranNH3 - r_nit + dNO3 <- tranNO3 + r_nit + dO2 <- tranO2 - 2 * r_nit + reaeration + + list(c( dNH3, dNO3, dO2 )) +} + +print(system.time( +std <- steady.1D(y = runif(3 * N), parms = NULL, + names=c("NH3", "NO3", "O2"), + func = Estuary, dimens = N, + positive = TRUE) +)) + +NH3in <- 100 +std2 <- steady.1D(y = runif(3 * N), parms = NULL, + names=c("NH3", "NO3", "O2"), + func = Estuary, dimens = N, + positive = TRUE) + +plot(std, std2, grid = Grid$x.mid, ylab = "mmol/m3", + xlab = "m", mfrow = c(1,3), col = "black") +legend("bottomright", lty = 1:2, title = "NH3in", + legend = c(500, 100)) diff --git a/Solving_Differential_Equations_in_R_Cash/DAEinR.R b/Solving_Differential_Equations_in_R_Cash/DAEinR.R new file mode 100644 index 000000000..f6bbca479 --- /dev/null +++ b/Solving_Differential_Equations_in_R_Cash/DAEinR.R @@ -0,0 +1,212 @@ +## ============================================================================= +## +## Examples from the book +## Soetaert, K., Cash, J.R. and Mazzia, F. (2012). +## Solving Differential Equations in R. +## Springer +## +## Chapter 6. Solving Differential Algebraic Equations in R. +## +## ============================================================================= + +## ----------------------------------------------------------------------------- +## A simple DAE +## ----------------------------------------------------------------------------- + +resdae <- function (t, y, dy, p) { + r1 <- dy[1] - y[2] + r2 <- y[1] - cos(t) + list(c(r1, r2)) +} + +library(deTestSet) +yini <- c(y1 = cos(0), y2 = -sin(0)) +dyini <- c(-sin(0), -cos(0)) +times <- seq(from = 0, to = 10, by = 0.1) +index <- c(1, 1, 0) +out1 <- mebdfi(times = times, res = resdae, y = yini, + atol = 1e-10, rtol = 1e-10, dy = dyini, + parms = NULL, nind = index) + +max (abs(out1[,"y1"] - cos(times)), abs(out1[,"y2"] + sin(times))) + +fundae <- function (t, y, p) { + f1 <- y[2] + f2 <- y[1] - cos(t) + list(c(f1, f2)) +} +M <- matrix(nrow = 2, ncol = 2, data = c(1, 0, 0, 0)) +out2 <- radau(times = times, fun = fundae, y = yini, + atol = 1e-10, rtol = 1e-10, mass = M, + parms = NULL, nind = index) +max (abs(out2[,"y1"] - cos(times)), abs(out2[,"y2"] + sin(times))) + + +## ----------------------------------------------------------------------------- +## implicit ODE +## ----------------------------------------------------------------------------- +implicit <- function(t, y, dy, parms) { + list(t*y^2*dy^3 - y^3*dy^2 + t*(t^2+1)*dy - t^2*y) +} +yini <- sqrt(3/2) +times <- seq(from = 1, to = 10, by = 0.1) + +library(rootSolve) +rootfun <- function (dy, y, t) + t*y^2*dy^3 - y^3*dy^2 + t*(t^2+1)*dy - t^2*y +dyini <- multiroot(f = rootfun, start = 0, y = yini, + t = times[1] )$root +dyini + +out <- mebdfi(times = times, res = implicit, y = yini, + dy = dyini, parms = NULL) +out2 <- daspk (times = times, res = implicit, y = yini, + dy = dyini, parms = NULL) + +max(abs(out [,2]- sqrt(times^2+0.5))) +max(abs(out2[,2]- sqrt(times^2+0.5))) + +implicit2 <- function (t, y, p) { + f1 <- y[2] + f2 <- t*y[1]^2*y[2]^3-y[1]^3*y[2]^2+t*(t^2+1)*y[2]-t^2*y[1] + list(c(f1, f2)) + } +M <- matrix(nrow = 2, ncol = 2, data = c(1, 0, 0, 0)) +yini_li <- c(yini,dyini) +out3 <- radau(times = times, fun = implicit2, y = yini_li, + mass = M, parms = NULL) +out4 <- gamd (times = times, fun = implicit2, y = yini_li, + mass = M, parms = NULL) +max(abs(out3[,2]- sqrt(times^2+0.5))) +max(abs(out4[,2]- sqrt(times^2+0.5))) + +## ----------------------------------------------------------------------------- +## The pendulum equation +## ----------------------------------------------------------------------------- +library(deTestSet) + +pendulum <- function (t, y, dy, parms) { + list(c(-dy[1] + y[3] , + -dy[2] + y[4] , + -dy[3] -y[5]*y[1] , + -dy[4] -y[5]*y[2] - 9.8, + y[1]^2 + y[2]^2 -1 + )) +} + +yini <- c(x = 1, y = 0, u = 0, v = 1 , lam = 1) +dyini <- c(dx = 0,dy = 1,du = -1,dv = -9.8,dlam = 3*9.8) +times <- seq(from = 0, to = 10, by = 0.01) +index3 <- c(2, 2, 1) +out3 <- mebdfi (y = yini, dy = dyini, res = pendulum, + parms = NULL, times = times, + nind = index3) + +plot(out3, lwd = 2) +plot(out3[, 2:3]) +mtext(side = 3, outer = TRUE, line = -1.5, + "Pendulum", cex = 1.5) + +## ----------------------------------------------------------------------------- +## The car axis problem +## ----------------------------------------------------------------------------- +caraxis <- function(t, y, dy, parms) { + with(as.list(y), { + f <- rep(0, 10) + yb <- r * sin(w * t) + xb <- sqrt(L^2 - yb^2) + Ll <- sqrt(xl^2 + yl^2) + Lr <- sqrt((xr - xb)^2 + (yr - yb)^2) + f[1:4] <- y[5:8] + f[5] <- 1/k*((L0-Ll)*xl/Ll + lam1*xb + 2*lam2*(xl-xr)) + f[6] <- 1/k*((L0-Ll)*yl/Ll + lam1*yb + 2*lam2*(yl-yr)) -g + f[7] <- 1/k*((L0-Lr)*(xr - xb)/Lr - 2*lam2*(xl-xr)) + f[8] <- 1/k*((L0-Lr)*(yr - yb)/Lr - 2*lam2*(yl-yr)) -g + f[9] <- xb * xl + yb * yl + f[10]<- (xl - xr)^2 + (yl - yr)^2 - L^2 + + delt <- dy - f + delt[9:10] <- -f[9:10] + + list(delt) + }) +} + +eps <- 0.01; M <- 10; k <- M * eps * eps/2 +L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 9.8 + +yini <- c(xl = 0, yl = L0, xr = L, yr = L0, + ul = -L0/L, vl = 0, ur = -L0/L, vr = 0, + lam1 = 0, lam2 = 0) + +library(rootSolve) +rootfun <- function (dyi, y, t) + unlist(caraxis(t, y, dy = c(dyi, 0, 0), + parms = NULL)) [1:8] + +dyini <- multiroot(f = rootfun, start = rep(0,8), + y = yini, t = 0)$root +(dyini <- c(dyini,0,0)) +caraxis(t = 0, yini, dyini, NULL) + +index <- c(4, 4, 2) +times <- seq(from = 0, to = 3, by = 0.01) +out <- mebdfi(y = yini, dy = dyini, times = times, + res = caraxis, parms = parameter, nind = index) + +par(mar = c(4, 4, 3, 2)) +plot(out, lwd = 2, mfrow = c(4,3)) +plot(out[,c("xl", "yl")], xlab = "xleft", ylab = "yleft", + type = "l", lwd = 2) +plot(out[,c("xr", "yr")], xlab = "xright", ylab = "yright", + type = "l", lwd = 2) + +## ----------------------------------------------------------------------------- +## The transistor equation +## ----------------------------------------------------------------------------- +library(deSolve) + +Transistor <- function(t, u, du, pars) { + delt <- vector(length = 8) + uin <- 0.1 * sin(200 * pi * t) + g23 <- beta * (exp( (u[2] - u[3]) / uf) - 1) + g56 <- beta * (exp( (u[5] - u[6]) / uf) - 1) + + delt[1] <- (u[1] - uin)/R0 + delt[2] <- u[2]/R1 + (u[2]-ub)/R2 + (1-alpha) * g23 + delt[3] <- u[3]/R3 - g23 + delt[4] <- (u[4] - ub) / R4 + alpha * g23 + delt[5] <- u[5]/R5 + (u[5]-ub)/R6 + (1-alpha) * g56 + delt[6] <- u[6]/R7 - g56 + delt[7] <- (u[7] - ub) / R8 + alpha * g56 + delt[8] <- u[8]/R9 + list(delt) +} + +ub <- 6; uf <- 0.026; alpha <- 0.99; beta <- 1e-6; R0 <- 1000 +R1 <- R2 <- R3 <- R4 <- R5 <- R6 <- R7 <- R8 <- R9 <- 9000 +C1 <- 1e-6; C2 <- 2e-6; C3 <- 3e-6; C4 <- 4e-6; C5 <- 5e-6 + +mass <- matrix(nrow = 8, ncol = 8, byrow = TRUE, data = c( + -C1,C1, 0, 0, 0, 0, 0, 0, + C1,-C1, 0, 0, 0, 0, 0, 0, + 0, 0,-C2, 0, 0, 0, 0, 0, + 0, 0, 0,-C3, C3, 0, 0, 0, + 0, 0, 0, C3,-C3, 0, 0, 0, + 0, 0, 0, 0, 0,-C4, 0, 0, + 0, 0, 0, 0, 0, 0,-C5, C5, + 0, 0, 0, 0, 0, 0, C5,-C5 +)) + +yini <- c(0, ub/(R2/R1+1), ub/(R2/R1+1), + ub, ub/(R6/R5+1), ub/(R6/R5+1), ub, 0) +names(yini) <- paste("u", 1:8, sep = "") + +ind <- c(8, 0, 0) +times <- seq(from = 0, to = 0.2, by = 0.001) + +out <- radau(func = Transistor, y = yini, parms = NULL, + times = times, mass = mass, nind = ind) + +plot(out, lwd = 2, which = c("u1", "u5", "u8"), + mfrow = c(1, 3)) diff --git a/Solving_Differential_Equations_in_R_Cash/DDEinR.R b/Solving_Differential_Equations_in_R_Cash/DDEinR.R new file mode 100644 index 000000000..51cc8d0db --- /dev/null +++ b/Solving_Differential_Equations_in_R_Cash/DDEinR.R @@ -0,0 +1,179 @@ +## ============================================================================= +## +## Examples from the book +## Soetaert, K., Cash, J.R. and Mazzia, F. (2012). +## Solving Differential Equations in R. +## Springer +## +## Chapter 8. Solving Delay Differential Equations in R. +## +## ============================================================================= + +## ----------------------------------------------------------------------------- +## Simple example 1 +## ----------------------------------------------------------------------------- +library(deSolve) +DDE1 <- function(t, y, parms) { + tlag <- t - 1 + if (tlag <= 0) + ylag <- 1 + else + ylag <- lagvalue(tlag) + + list(dy = - ylag, ylag = ylag) +} +yinit <- 1 +times <- seq(from = 0, to = 10, by = 0.1) +yout <- dede(y = yinit, times = times, func = DDE1, + parms = NULL, atol = 1e-10, rtol = 1e-10 ) + +tt <- which(times >= 1 & times <= 2) +analytic <- c(1-times[times <1] , 0.5*times[tt]^2 - 2*times[tt]+3/2) +max(abs(yout[times <= 2,2] - analytic)) + +## ----------------------------------------------------------------------------- +## Simple example 2 +## ----------------------------------------------------------------------------- +DDE2 <- function(t, y, parms) { + tlag <- t - 1 + if (tlag <= 0) + ylag <- 1 + else + ylag <- lagderiv(tlag) + + list(dy = - ylag, ylag = ylag) +} +yout2 <- dede(y = yinit, times = times, func = DDE2, + parms = NULL ) + + +## ----------------------------------------------------------------------------- +## White Blood cells +## ----------------------------------------------------------------------------- +library(deSolve) + +mackey <- function(t, y, parms, tau) { + tlag <- t - tau + if (tlag <= 0) + ylag <- 0.5 + else + ylag <- lagvalue(tlag) + dy <- 0.2 * ylag * 1/(1+ylag^10) - 0.1 * y + list(dy = dy, ylag = ylag) +} + +yinit <- 0.5 +times <- seq(from = 0, to = 300, by = 0.1) + +yout1 <- dede(y = yinit, times = times, func = mackey, + parms = NULL, tau = 10) +yout2 <- dede(y = yinit, times = times, func = mackey, + parms = NULL, tau = 20) + +plot(yout1, lwd = 2, main = "tau=10", + ylab = "y", mfrow = c(2, 2), which = 1) +plot(yout1[,-1], type = "l", lwd = 2, xlab = "y") +plot(yout2, lwd = 2, main = "tau=20", + ylab = "y", mfrow = NULL, which = 1) +plot(yout2[,-1], type = "l", lwd = 2, xlab = "y") + +## ----------------------------------------------------------------------------- +## DDE with root +## ----------------------------------------------------------------------------- +xb <- -0.427; a <- 0.16; xi <- 0.02; u <- 0.5; tau <- 1 +yinit <- c(y = 0.6) + +mariott <- function(t, y, parms) { + tlag <- t - 12 + if (tlag <= 0) + ylag <- 0.6 + else + ylag <- lagvalue(tlag) + + Delt <- ylag - xb + sDelt <- sign(Delt) + + dy <- (-y + pi*(a + xi*sDelt - u*(sin(Delt))^2))/tau + list(dy) +} + +times <- seq(from = 0, to = 120, by = 0.5) +yout <- dede(y = yinit, times = times, func = mariott, + parms = NULL) + +root <- function(t, y, parms) { + tlag <- t - 12 + if (tlag <= 0) + return (1) # not a root + else + return(lagvalue(tlag)- xb) +} + +event <- function(t, y, parms) return(y) + +yout <- dede(y = yinit, times = times, func = mariott, + parms = NULL, rootfun = root, + events = list(func = event, root = TRUE)) + +attributes(yout)$troot + +plot(yout, lwd = 2, + main = "Controller problem") +abline(v = attributes(yout)$troot, col = "grey") + +## ----------------------------------------------------------------------------- +## Vanishing time delay +## ----------------------------------------------------------------------------- +vanishing <- function(t, y, parms, cc) { + + tlag <- t*y^2 + if (tlag <= 0) { + ylag <- 0 + dylag <- 0 + } else { + ylag <- lagvalue(tlag) + dylag <- lagderiv(tlag) + } + dy <- cos(t)*(1+ylag) + cc*y*dylag + + (1-cc)*sin(t)*cos(t*sin(t)^2) - sin(t+t*sin(t)^2) + + list(dy) +} +yinit <- c(y = 0) +times <- seq(from = 0, to = 2*pi, by = 0.1) +yout <- dede(y = 0, times = times, func = vanishing, + parms = NULL, cc = -0.5, + atol = 1e-10, rtol = 1e-10) +print(max(abs(yout[,2] - sin(yout[,1])))) + +## ----------------------------------------------------------------------------- +## Predator-prey with harvesting +## ----------------------------------------------------------------------------- + +LVdede <- function(t, y, p) { + if (t > tau1) Lag1 <- lagvalue(t - tau1) else Lag1 <- yini + if (t > tau2) Lag2 <- lagvalue(t - tau2) else Lag2 <- yini + + dy1 <- r * y[1] *(1 - Lag1[1]/K) - a*y[1]*y[2] + dy2 <- a * b * Lag2[1]*Lag2[2] - d*y[2] + + list(c(dy1, dy2)) +} + +rootfun <- function(t, y, p) + return(y[1] - Ycrit) + +eventfun <- function(t, y, p) + return (c(y[1] * 0.7, y[2])) + +r <- 1; K <- 1; a <- 2; b <- 1; d <- 1; Ycrit <- 1.2*d/(a*b) +tau1 <- 0.2; tau2 <- 0.2 + +yini <- c(y1 = 0.2, y2 = 0.1) +times <- seq(from = 0, to = 200, by = 0.01) +yout <- dede(func = LVdede, y = yini, times = times, + parms = 0, rootfun = rootfun, + events = list(func = eventfun, root = TRUE)) +attributes(yout)$troot [1:10] + +plot(yout[,-1], type = "l") diff --git a/Solving_Differential_Equations_in_R_Cash/ODEinR.R b/Solving_Differential_Equations_in_R_Cash/ODEinR.R new file mode 100644 index 000000000..b75750bc1 --- /dev/null +++ b/Solving_Differential_Equations_in_R_Cash/ODEinR.R @@ -0,0 +1,340 @@ +## ============================================================================= +## +## Examples from the book +## Soetaert, K., Cash, J.R. and Mazzia, F. (2012). +## Solving Differential Equations in R. +## Springer +## +## Chapter 4. Solving Ordinary Differential Equations in R. +## +## ============================================================================= + + +## ----------------------------------------------------------------------------- +## One variable +## ----------------------------------------------------------------------------- +r <- 1 ; K <- 10 +yini <- c(y = 2) + +derivs <- function(t, y, parms) + list(r * y * (1-y/K)) + +library(deSolve) +times <- seq(from = 0, to = 20, by = 0.2) +out <- ode(y = yini, times = times, func = derivs, + parms = NULL) + +head(out, n = 3) + +# +yini <- c(y = 12) +out2 <- ode(y = yini, times = times, func = derivs, + parms = NULL) + +plot(out, out2, main = "logistic growth", lwd = 2) + +## ----------------------------------------------------------------------------- +## The Lorenz equations +## ----------------------------------------------------------------------------- +a <- -8/3 ; b <- -10; c <- 28 +yini <- c(X = 1, Y = 1, Z = 1) + +Lorenz <- function (t, y, parms) { + with(as.list(y), { + dX <- a * X + Y * Z + dY <- b * (Y - Z) + dZ <- -X * Y + c * Y - Z + list(c(dX, dY, dZ)) + }) +} + +times <- seq(from = 0, to = 100, by = 0.01) +out <- ode(y = yini, times = times, func = Lorenz, parms = NULL) + +plot(out, lwd = 2) +plot(out[,"X"], out[,"Y"], type = "l", xlab = "X", + ylab = "Y", main = "butterfly") + +## ----------------------------------------------------------------------------- +## The rigid ODE +## ----------------------------------------------------------------------------- +library(deSolve) +yini <- c(1, 0, 0.9) + +rigidode <- function(t, y, parms) { + dy1 <- -2 * y[2] * y[3] + dy2 <- 1.25* y[1] * y[3] + dy3 <- -0.5* y[1] * y[2] + list(c(dy1, dy2, dy3)) +} + +times <- seq(from = 0, to = 20, by = 0.01) +out <- ode (times = times, y = yini, func = rigidode, + parms = NULL, method = rkMethod("rk45ck")) +head (out, n = 3) + +par(mfrow = c(1, 2)) +matplot(x = out[,1], y = out[,-1], type = "l", lwd = 2, + lty = "solid", col = c("red", "blue", "black"), + xlab = "times", ylab = "y", main = "rigidode") + +legend("bottomright", col = c("red", "blue", "black"), + legend = c("y1", "y2", "y3"), lwd = 2) + +library(scatterplot3d) +par(mar = c(0,0,0,0)) +scatterplot3d(out[,-1], type = "l", lwd = 2, xlab = "", + ylab = "", zlab = "", main = "rigidode") + +## ----------------------------------------------------------------------------- +## Arenstorff orbits +## ----------------------------------------------------------------------------- +library(deSolve) +Arenstorff <- function(t, y, p) { + D1 <- ((y[1] + mu1)^2 + y[2]^2)^(3/2) + D2 <- ((y[1] - mu2)^2 + y[2]^2)^(3/2) + dy1 <- y[3] + dy2 <- y[4] + dy3 <- y[1] + 2*y[4] - mu2*(y[1]+mu1)/D1 - mu1*(y[1]-mu2)/D2 + dy4 <- y[2] - 2*y[3] - mu2*y[2]/D1 - mu1*y[2]/D2 + return(list( c(dy1, dy2, dy3, dy4) )) +} +mu1 <- 0.012277471 +mu2 <- 1 - mu1 +yini <- c(y1 = 0.994, y2 = 0, + dy1 = 0, dy2 = -2.00158510637908252240537862224) +times <- seq(from = 0, to = 18, by = 0.01) + +out <- ode(func = Arenstorff, y = yini, times = times, + parms = 0, method = "ode45") + +yini2 <- c(y1 = 0.994, y2 = 0, + dy1 = 0, dy2 = -2.0317326295573368357302057924) +out2 <- ode(func = Arenstorff, y = yini2, times = times, + parms = 0, method = "ode45") + +yini3 <- c(y1 = 1.2, y2 = 0, + dy1 = 0, dy2 = -1.049357510) +out3 <- ode(func = Arenstorff, y = yini3, times = times, + parms = 0, method = "ode45") + +plot(out, out2, out3, which = c("y1", "y2"), + mfrow = c(2, 2), col = "black", lwd = 2) + +plot(out[ ,c("y1", "y2")], type = "l", lwd = 2, + xlab = "y1", ylab = "y2", main = "solutions 1,2") +lines(out2[ ,c("y1", "y2")], lwd = 2, lty = 2) + +plot(out3[ ,c("y1", "y2")], type = "l", lwd = 2, lty = 3, + xlab = "y1", ylab = "y2", main = "solution 3") + +## ----------------------------------------------------------------------------- +## The pleiades +## ----------------------------------------------------------------------------- +library(deSolve) + +pleiade <- function (t, Y, pars) { + x <- Y[1:7] + y <- Y[8:14] + u <- Y[15:21] + v <- Y[22:28] + + distx <- outer(x, x, FUN = function(x, y) x - y) + disty <- outer(y, y, FUN = function(x, y) x - y) + + rij3 <- (distx^2 + disty^2)^(3/2) + + fx <- starMass * distx / rij3 + fy <- starMass * disty / rij3 + + list(c(dx = u, + dy = v, + du = colSums(fx, na.rm = TRUE), + dv = colSums(fy, na.rm = TRUE))) +} + +starMass <- 1:7 + +yini<- c(x1= 3, x2= 3, x3=-1, x4=-3, x5= 2, x6=-2, x7= 2, + y1= 3, y2=-3, y3= 2, y4= 0, y5= 0, y6=-4, y7= 4, + u1= 0, u2= 0, u3= 0, u4= 0, u5= 0, u6=1.75, u7=-1.5, + v1= 0, v2= 0, v3= 0, v4=-1.25, v5= 1, v6= 0, v7= 0) +print(system.time( +out <- ode(func = pleiade, parms = NULL, y = yini, + method = "adams", times = seq(0, 3, 0.01)))) + +par(mfrow = c(3, 3)) +for (i in 1:7) { + plot(out[,i+1], out[,i+8], type = "l", + main = paste("star ",i), xlab = "x", ylab = "y") + points (yini[i], yini[i+7]) +} +plot(out[, 2:8], out[, 9:15], type = "p", cex = 0.5, + main = "ALL", xlab = "x", ylab = "y") +text(yini[1:7], yini[8:14], 1:7) +matplot(out[,"time"], out[, c("u1", "u7")], type = "l", + lwd = 2, col = c("black", "grey"), lty = 1, + xlab = "time", ylab = "velocity", main = "stars 1, 7") +abline(v = c(1.23, 1.68), lty = 2) +legend("bottomright", col = c("black", "grey"), lwd = 2, + legend = c("u1", "u7")) + +## ----------------------------------------------------------------------------- +## The ozone model +## ----------------------------------------------------------------------------- +load(file = "Light.rda") +head(Light, n = 4) + +irradiance <- approxfun(Light) +irradiance(seq(from = 0, to = 1, by = 0.25)) +k3 <- 1e-11; k2 <- 1e10; k1a <- 1e-30 +k1b <- 1; sigma <- 1e11 +yini <- c(O = 0, NO = 1.3e8, NO2 = 5e11, O3 = 8e11) + +chemistry <- function(t, y, parms) { + with(as.list(y), { + + radiation <- irradiance(t) + k1 <- k1a + k1b*radiation + + dO <- k1*NO2 - k2*O + dNO <- k1*NO2 - k3*NO*O3 + sigma + dNO2 <- -k1*NO2 + k3*NO*O3 + dO3 <- k2*O - k3*NO*O3 + list(c(dO, dNO, dNO2, dO3), radiation = radiation) + }) +} + +times <- seq(from = 8/24, to = 5, by = 0.01) +out <- ode(func = chemistry, parms = NULL, y = yini, + times = times, method = "bdf") + +plot(out, type = "l", lwd = 2 ) + +## ----------------------------------------------------------------------------- +## Pharmacokinetic 1 +## ----------------------------------------------------------------------------- + +a <- 6; b <- 0.6 +yini <- c(intestine = 0, blood = 0) + +pharmacokinetics <- function(t, y, p) { + if ( (24*t) %% 24 <= 1) + uptake <- 2 + else + uptake <- 0 + dy1 <- - a* y[1] + uptake + dy2 <- a* y[1] - b *y[2] + list(c(dy1, dy2)) +} + +times <- seq(from = 0, to = 10, by = 1/24) +out <- ode(func = pharmacokinetics, times = times, + y = yini, parms = NULL) + +times <- seq(0, 10, by = 3/24) +out2 <- ode(func = pharmacokinetics, times = times, + y = yini, parms = NULL, method = "impAdams") + +plot(out, lwd = 2, mfrow = c(2, 2), xlab = "day") +plot(out2, lwd = 2, mfrow = NULL, xlab = "day") + +## ----------------------------------------------------------------------------- +## Pharmacokinetic 2 +## ----------------------------------------------------------------------------- +b <- 0.6 +yini <- c(blood = 0) + +pharmaco2 <- function(t, blood, p) { + dblood <- - b * blood + list(dblood) +} + +injectevents <- data.frame(var = "blood", + time = 0:20, + value = 40, + method = "add") +head(injectevents, n=3) + +times <- seq(from = 0, to = 10, by = 1/24) + +out2 <- ode(func = pharmaco2, times = times, y = yini, + parms = NULL, method = "impAdams", + events = list(data = injectevents)) + +plot(out2, lwd = 2) + +## ----------------------------------------------------------------------------- +## Bouncing ball +## ----------------------------------------------------------------------------- +library(deSolve) +yini <- c(height = 0, velocity = 10) + +ball <- function(t, y, parms) { + dy1 <- y[2] + dy2 <- -9.8 + + list(c(dy1, dy2)) +} + +rootfunc <- function(t, y, parms) y[1] + +eventfunc <- function(t, y, parms) { + y[1] <- 0 + y[2] <- -0.9*y[2] + return(y) +} + +times <- seq(from = 0, to = 20, by = 0.01) +out <- ode(times = times, y = yini, func = ball, + parms = NULL, rootfun = rootfunc, + events = list(func = eventfunc, root = TRUE)) + +plot(out, which = "height", lwd = 2, + main = "bouncing ball", ylab = "height") + +## ----------------------------------------------------------------------------- +## Temperature in a room +## ----------------------------------------------------------------------------- +yini <- c(temp = 18, heating_on = 1) + +temp <- function(t, y, parms) { + dy1 <- ifelse(y[2] == 1, 1.0, -0.5) + dy2 <- 0 + list(c(dy1, dy2)) +} + +rootfunc <- function(t, y, parms) c(y[1]-18, y[1]-20) + +eventfunc <- function(t, y, parms) { + y[1] <- y[1] + y[2] <- ! y[2] + return(y) +} + +times <- seq(from = 0, to = 20, by = 0.1) + +out <- lsode(times = times, y = yini, func = temp, + parms = NULL, rootfun = rootfunc, + events = list(func = eventfunc, root = TRUE)) + +attributes(out)$troot +plot(out, lwd = 2) + +## ----------------------------------------------------------------------------- +## The van der Pol equation +## ----------------------------------------------------------------------------- + +yini <- c(y = 2, dy = 0) + +Vdpol <- function(t, y, mu) + list(c(y[2], + mu * (1 - y[1]^2) * y[2] - y[1])) + +times <- seq(from = 0, to = 30, by = 0.01) + +nonstiff <- ode(func = Vdpol, parms = 1, y = yini, + times = times) +diagnostics(nonstiff) + + diff --git a/Solving_Differential_Equations_in_R_Cash/PDEinR.R b/Solving_Differential_Equations_in_R_Cash/PDEinR.R new file mode 100644 index 000000000..e4e5f8145 --- /dev/null +++ b/Solving_Differential_Equations_in_R_Cash/PDEinR.R @@ -0,0 +1,352 @@ +## ============================================================================= +## +## Examples from the book +## Soetaert, K., Cash, J.R. and Mazzia, F. (2012). +## Solving Differential Equations in R. +## Springer +## +## Chapter 10. Solving Partial Differential Equations in R. +## +## ============================================================================= + +## ----------------------------------------------------------------------------- +## The heat equation +## ----------------------------------------------------------------------------- +library(ReacTran) +N <- 100 +xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) +x <- xgrid$x.mid +D.coeff <- 0.01 + +Diffusion <- function (t, Y, parms){ + tran <- tran.1D(C = Y, C.up = 0, C.down = 1, + D = D.coeff, dx = xgrid) + list(dY = tran$dC, flux.up = tran$flux.up, + flux.down = tran$flux.down) +} + +Yini <- sin(pi*x) +times <- seq(from = 0, to = 5, by = 0.01) +print(system.time( + out <- ode.1D(y = Yini, times = times, func = Diffusion, + parms = NULL, dimens = N) )) + +par (mfrow=c(1, 2)) +plot(out[1, 2:(N+1)], x, type = "l", lwd = 2, + xlab = "Variable, Y", ylab = "Distance, x") +for (i in seq(2, length(times), by = 50)) + lines(out[i, 2:(N+1)], x) +image(out, grid = x, mfrow = NULL, ylab = "Distance, x", + main = "Y") + +## ----------------------------------------------------------------------------- +## The wave equation +## ----------------------------------------------------------------------------- +library(ReacTran) +dx <- 0.2 +xgrid <- setup.grid.1D(x.up = -100, x.down = 100, dx.1 = dx) +x <- xgrid$x.mid +N <- xgrid$N + +lam <- 0.05 +uini <- exp(-lam*x^2) +vini <- rep(0, N) +yini <- c(uini, vini) +times <- seq (from = 0, to = 50, by = 1) + +wave <- function (t, y, parms) { + u <- y[1:N] + v <- y[(N+1):(2*N)] + + du <- v + dv <- tran.1D(C = u, C.up = 0, C.down = 0, D = 1, + dx = xgrid)$dC + + return(list(c(du, dv))) +} + +out <- ode.1D(func = wave, y = yini, times = times, + parms = NULL, method = "adams", + dimens = N, names = c("u", "v")) +u <- subset(out, which = "u") + +analytic <- function (t, x) + 0.5 * (exp(-lam * (x+1*t)^2 ) +exp(-lam * (x-1*t)^2) ) + +OutAna <- outer(times, x, FUN = analytic) + +max(abs(u - OutAna)) + +outtime <- seq(from = 0, to = 50, by = 10) +matplot.1D(out, which = "u", subset = time %in% outtime, + grid = x, xlab = "x", ylab = "u", type = "l", + lwd = 2, xlim = c(-50, 50), + col = c("black", rep("darkgrey", 5))) +legend("topright", lty = 1:6, lwd = 2, + col = c("black", rep("darkgrey", 5)), + title = "t = ", legend = outtime) + +## ----------------------------------------------------------------------------- +## The Laplace equation +## ----------------------------------------------------------------------------- +Nx <- 100 +Ny <- 100 +xgrid <- setup.grid.1D (x.up = 0, x.down = 1, N = Nx) +ygrid <- setup.grid.1D (x.up = 0, x.down = 1, N = Ny) +x <- xgrid$x.mid +y <- ygrid$x.mid + +laplace <- function(t, U, parms) { + w <- matrix(nrow = Nx, ncol = Ny, data = U) + dw <- tran.2D(C = w, C.x.up = 0, C.x.down = 0, + flux.y.up = 0, + flux.y.down = -1 * sin(pi*x)*pi*sinh(pi), + D.x = 1, D.y = 1, + dx = xgrid, dy = ygrid)$dC + list(dw) +} + +print(system.time( + out <- steady.2D(y = runif(Nx*Ny), func = laplace, + parms = NULL, nspec = 1, + dimens = c(Nx, Ny), lrw = 1e7) +)) +w <- matrix(nrow = Nx, ncol = Ny, data = out$y) + +analytic <- function (x, y) sin(pi*x) * cosh(pi*y) +OutAna <- outer(x, y, FUN = analytic) + +max(abs(w - OutAna)) + +image(out, grid = list(x, y), main = "elliptic Laplace", + add.contour = TRUE) + +## ----------------------------------------------------------------------------- +## The advection equation +## ----------------------------------------------------------------------------- +adv.func <- function(t, y, p, adv.method) + list(advection.1D(C = y, C.up = y[N], C.down = y[1], + v = 0.1, adv.method = adv.method, + dx = xgrid)$dC) + +xgrid <- setup.grid.1D(0.3, 1.3, N = 50) +x <- xgrid$x.mid +N <- length(x) + +yini <- sin(pi * x)^50 +times <- seq(0, 20, 0.01) + +out1 <- ode.1D(y = yini, func = adv.func, times = times, + parms = NULL, method = "euler", dimens = N, + adv.method = "muscl") +out2 <- ode.1D(y = yini, func = adv.func, times = times, + parms = NULL, method = "euler", dimens = N, + adv.method = "super") +## ----------------------------------------------------------------------------- +## The Brusselator in 1D +## ----------------------------------------------------------------------------- + +library(ReacTran) +N <- 50 +Grid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) + +x1ini <- 1 + sin(2 * pi * Grid$x.mid) +x2ini <- rep(x = 3, times = N) +yini <- c(x1ini, x2ini) + +brusselator1D <- function(t, y, parms) { + + X1 <- y[1:N] + X2 <- y[(N+1):(2*N)] + + dX1 <- 1 + X1^2*X2 - 4*X1 + + tran.1D (C = X1, C.up = 1, C.down = 1, + D = 0.02, dx = Grid)$dC + dX2 <- 3*X1 - X1^2*X2 + + tran.1D (C = X2, C.up = 3, C.down = 3, + D = 0.02, dx = Grid)$dC + + list(c(dX1, dX2)) +} + +times <- seq(from = 0, to = 10, by = 0.1) + +print(system.time( + out <- ode.1D(y = yini, func = brusselator1D, + times = times, parms = NULL, nspec = 2, + names = c("X1", "X2"), dimens = N) +)) + +par(mfrow = c(2, 2)) +image(out, mfrow = NULL, grid = Grid$x.mid, + which = "X1", method = "contour") +image(out, mfrow = NULL, grid = Grid$x.mid, + which = "X1") + +par(mar = c(1, 1, 1, 1)) +image(out, mfrow = NULL, grid = Grid$x.mid, + which = "X1", method = "persp", col = NA) +image(out, mfrow = NULL, grid = Grid$x.mid, + which = "X1", method = "persp", border = NA, + shade = 0.3 ) + +## ----------------------------------------------------------------------------- +## The Brusselator in 2D +## ----------------------------------------------------------------------------- +brusselator2D <- function(t, y, parms) { + + X1 <- matrix(nrow = Nx, ncol = Ny, + data = y[1:(Nx*Ny)]) + X2 <- matrix(nrow = Nx, ncol = Ny, + data = y[(Nx*Ny+1) : (2*Nx*Ny)]) + + dX1 <- 1 + X1^2*X2 - 4*X1 + + tran.2D (C = X1, D.x = D_X1, D.y = D_X1, + dx = Gridx, dy = Gridy)$dC + dX2 <- 3*X1 - X1^2*X2 + + tran.2D (C = X2, D.x = D_X2, D.y = D_X2, + dx = Gridx, dy = Gridy)$dC + + list(c(dX1, dX2)) +} + +library(ReacTran) +Nx <- 50 +Ny <- 50 +Gridx <- setup.grid.1D(x.up = 0, x.down = 1, N = Nx) +Gridy <- setup.grid.1D(x.up = 0, x.down = 1, N = Ny) + +D_X1 <- 2 +D_X2 <- 8*D_X1 + +X1ini <- matrix(nrow = Nx, ncol = Ny, data = runif(Nx*Ny)) +X2ini <- matrix(nrow = Nx, ncol = Ny, data = runif(Nx*Ny)) + +yini <- c(X1ini, X2ini) + +times <- 0:8 +print(system.time( + out <- ode.2D(y = yini, parms = NULL, func = brusselator2D, + nspec = 2, dimens = c(Nx, Ny), times = times, + lrw = 2000000, names=c("X1", "X2")) +)) + +par(oma = c(0,0,1,0)) +image(out, which = "X1", xlab = "x", ylab = "y", + mfrow = c(3, 3), ask = FALSE, + main = paste("t = ", times), + grid = list(x = Gridx$x.mid, y = Gridy$x.mid)) +mtext(side = 3, outer = TRUE, cex = 1.25, line = -1, + "2-D Brusselator, species X1") + +## ----------------------------------------------------------------------------- +## The polar Laplace equation +## ----------------------------------------------------------------------------- +library(ReacTran) + +Nr <- 100 +Np <- 100 +r <- seq(2, 4, len = Nr+1) +theta <- seq(0, 2*pi, len = Np+1) +theta.mid <- 0.5*(theta[-1] + theta[-Np]) + +Model <- function(t, C, p) { + y = matrix(nrow = Nr, ncol = Np, data = C) + tran <- tran.polar (y, D.r = 1, r = r, theta = theta, + C.r.up = 0, C.r.down = 4 * sin(5*theta.mid), + cyclicBnd = 2) + list(tran$dC) +} + +STD <- steady.2D(y = runif(Nr*Np), parms = NULL, + func = Model, dimens = c(Nr, Np), + lrw = 1e6, cyclicBnd = 2) + +OUT <- polar2cart (STD, r = r, theta = theta, + x = seq(-4, 4, len = 400), + y = seq(-4, 4, len = 400)) + +image(OUT, main = "Laplace") + +## ----------------------------------------------------------------------------- +## The sine-gordon +## ----------------------------------------------------------------------------- + +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) + +## ----------------------------------------------------------------------------- +## The Schrodinger equation +## ----------------------------------------------------------------------------- +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") diff --git a/Solving_Differential_Equations_in_R_Cash/light.rda b/Solving_Differential_Equations_in_R_Cash/light.rda new file mode 100644 index 0000000000000000000000000000000000000000..f76111d9d9ad024c6dc6bec94930eee7b53626cc GIT binary patch literal 1719 zcmZ|DeKgYx0|#&_hRu4)ZPGH=BQI-SZY?(~bgpo{gvD7%{ZzVSF_O1z6t11TgiAc@ zCPZEm>mt$^JGFa-jIDV~+awL6%}ks3$Mr|M|J?8AobUObuR)fE%67oq;;VjsR}PY` z(0wJOjm;_KmyHmID3yZV$qWC zYbglbD38;-f}-HKJr}D`Yl=KzT`*%W?C&Fei77=dK9u+OEtPy08!r|gffxnvS}Qo& z-FX8h+3mGvGIe~g?^)i!K8Vq$n};&63+K>j6TXwfXk!2Tw4Q@*msKjTqZ5*>wsu%y zEy0ganF}JpFhkTOVu%>8ii=r&hDZkNIB_|(5Ss;59U z_P2w#!P~ctLnA@0X_zj3!QBg9{^ZvQvOe&Uj=%2Q0r6FNG@mz`bRD9p!b?&jz%X>X zzL6XG7Y7M>~9Zd{X^h-)6Ft2nON}&@?k^qG{5#|5}ApfufP|^aF)BzTJE zY&p|`*d8D4xR0pL!26Nj*lE$gaN?~BsuMt_LrMB8cCziOob@lLu0SpwL9$=DDO;Kr z4Gvd)J1nPZW#Uvic%@?`H*q$;jo=eJJEOuS#~`oy;)nV#+cpn55XLpW8Zxa5izDJ} z!Xkc4jK}^jhwVBRcL!Z1Wu0@Luy#T&ujN%a8Lj3o&$*cSR;# z?$p*AL4{dlVFwNmFND^=zho9SQ27Hj+7#LFPR8_f1v9x6zqTKd%uZ`DQyza0dnR$L ztN_%4Wk(;+%}Ls8HXrh}h_I64VM3XVh~YvT3l(v&4SN3d+*4wG9hT1}eEM65UmuT0 zOgYkOalYZP`SNM)-Q9i{AIR44X!ZEVdthslZhM|@0cdlEA%SfJb2MuHUi?{;l-?S% z6^JMh@iIHr33PXJEXX~JUH5b-<5=3}$$3dURG-9Kls`?&#!4N1zoN`OvIjp#>m!=& zR`1Vz{A8JVj~d&GxmFl%b2zjbQLV@|;sk=_XV|-#wcf7lKO;BOe@C^1y=i`S>{CxA z&OLdm>{5L1l^*^mA!&T=P|pI9>QE#Ojr(kVpK?+iactqAC5M?ax`zJXvCr{icJ zc-~j5UtUP&eqa8MQiYCAVLev7VdSfVgxE%5I@>E_QvO%-8?tBRR2Dr5eVM=KD9`G; z?(>r4hfXB4+RcpD2<=l?cHJ+6?}{O!-QU&sHeiS$wib0K2#L)zs7#C^De1)|lyXK( ziWo~hO)w=Rt+Co3icM#X!T1e zSlp|5#_uBqg&T~AKYL%)H12n!FM2F^H!}?F_Yh+YriAaH{Q{w5=<(CVOTCjRR=Cr) zFu%)k24BoEzfdATvYMEUi3zV>R@wtz2``W$^x< z@ltsf(#F^?sP4ka+tsdL41Q#AgKOiy1P+i^6IuMG+)bXR>j9(GN6Wh28!ehDzux|n zwR)HA?5wx5)Xwd|^SP!$lBvV;=4gy~BFe_kdsvqYuUMTW`M|=0E*{iB$Sh@v?b|vU PH^K1Luop|aRaE{DvtV77 literal 0 HcmV?d00001