0

我正在尝试用一个 ODE 和一个二阶 PDE 求解一个微分方程组。我已经离散化了 PDE，以便我可以使用 deSolve ode.1D 来解决它，但现在我无法弄清楚如何将我的边界条件包含在空间 dE/dx 导数中。在 x = 0 和 x = 1 处，dE/dx 的边界值对于获得正确解很重要。如何将它们包含在模型中？

我的代码：

BVmod1D <- function(time, state, parms, N, dx) {
  with(as.list(parms), {
    E <- state[1 : N]
    U <- state[(N + 1) : (2 * N)]
    

    coef1 <- Sv * io / (Qox - Qred)
    coef2 <- Tau * io / Cd
    BV <- exp(aa * f * (E - 0.5 * (1 + 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 # dU/dt, Eqn 8
    dE <- (ddEddx - (coef2 * BV)) / Tau #  dE/dt, Eqn 12

    return(list(c(dE, dU)))
  })
}


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)] 
times <- seq(0, 1 , by = 0.002 ) 


tail(out <- ode.1D(
  y = yini,
  times = times,
  func = BVmod1D,
  parms = pars,
  nspec = 2,
  N = N,
  dx = dx
))
4

0 回答 0