R help - Apr 2021 - 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
)
```

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 at 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 at 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]]