How to force boundary conditions on discretized derivative in deSolve?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

How to force boundary conditions on discretized derivative in deSolve?

Brockway, Molly
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.
Reply | Threaded
Open this post in threaded view
|

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

Bert Gunter-2
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 <[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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.