[R] How to force boundary conditions on discretized derivative in deSolve?

[Bert Gunter]

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 <mbrockway using 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 using 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]]