[R] 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
)
```