# How to force boundary conditions on discretized derivative in deSolve?

2 messages
Open this post in threaded view
|

## How to force boundary conditions on discretized derivative in deSolve?

 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 ) ``` ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: How to force boundary conditions on discretized derivative in deSolve?

 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 ), 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 <[hidden email]> 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 > ) > ``` > > > > > > > > > > > > > > > > > > > ______________________________________________ > [hidden email] 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]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.