1. Your query is off topic. See the posting guide linked below for what is appropriate. In particular, note: "*Questions about statistics:* The R mailing lists are primarily intended for questions and discussion about the R software. However, questions about statistical methodology are sometimes posted. If the question is well-asked and of interest to someone on the list, it *may* elicit an informative up-to-date answer."
and "For questions about functions in standard packages distributed with R (see the FAQ Add-on packages in R <https://cran.r-project.org/doc/FAQ/R-FAQ.html#Add-on-packages-in-R>), ask questions on R-help. If the question relates to a *contributed package* , e.g., one downloaded from CRAN, try contacting the package maintainer first. You can also use find("functionname") and packageDescription("packagename") to find this information. *Only* send such questions to R-help or R-devel if you get no reply or need further assistance. This applies to both requests for help and to bug reports." 2. See https://cran.r-project.org/web/views/DifferentialEquations.html . Note especially the link to the dynamic models SIG therein, which I assume might be helpful to you. Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Wed, Apr 14, 2021 at 7:05 AM Brockway, Molly <mbrock...@mtech.edu> wrote: > Hello, > > I am using R package deSolve to solve a system of two differential > equations for a one-dimensional spatial and time-based problem. There is > one ODE and a second-order PDE. In order to solve with the function ode.1D, > I've discretized the spatial derivative and put both equations in terms of > the time derivative only. However, I'm now stuck with the problem of being > unable to impose boundary conditions on the spatial derivative for flux at > the edges of the system. How can I force a value for the discretized dE/dx > part of my equations at x = 0 and x = 1? > > The code I have been using is below. The ode.1D solver will run on it, but > the solutions aren't correct for my system owing to the missing boundary > conditions. > > Thanks, > > Molly C Brockway > Materials Science - PhD Candidate > Metallurgical and Materials Engineering - B.S. > > Montana Technological University > > > ``` > BVmod1D <- function(time, state, parms, N, dx) { > with(as.list(parms), { > U <- state[1 : N] > E <- state[(N + 1) : (2 * N)] > > coef1 <- Sv * io / (Qox - Qred) > coef2 <- Tau * io / Cd > BV <- exp(aa * f * (E - 0.5 * (U) )) - exp(-ac * f * (E - 0.5 * (1 + > U))) > > dEdx <- diff(c(E, E[N])) / dx # first spatial derivative of E > ddEddx <- diff(c(dEdx, dEdx[N])) / dx # second spatial derivative of E > > dU <- coef1 * BV > dE <- (ddEddx - (coef2 * BV)) / Tau > > return(list(c(dU, dE))) > }) > } > > pars <- c(Sv = 1.72e5, > io = 1.71e-6, > Qox = 1.56, > Qred = 0, > Tau = 3.79e-6, > Cd = 7.42e-8, > aa = 0.7674, > ac = 0.2326, > ks = 0.042992, > sig = 1.67e-6, > f = 38.92237) > > R <- 1 > N <- 50 > dx <- R / N > Vo <- 0.5 > > # Initial and Boundary Conditions > yini <- rep(0, (2 * N)) > yini[ 1 : N ] <- 2 * Vo > yini[ (N + 1) : (2 * N)] <- 1 > times <- seq(0, 1 , by = 0.002 ) > > out <- ode.1D( > y = yini, > times = times, > func = BVmod1D, > parms = pars, > nspec = 2, > N = N, > dx = dx > ) > ``` > > > > > > > > > > > > > > > > > > > ______________________________________________ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
