I'm having trouble figuring out how to implement boundary conditions for this problem: \begin{align} \frac{\partial n}{\partial t} &= D_n\nabla^2n - \nabla\cdot\left(\frac{\chi}{1+\alpha c}n\nabla c\right) - \nabla\cdot (pn\nabla f) \\ \frac{\partial f}{\partial t} &= \beta n-\gamma nf \\ \frac{\partial c}{\partial t} &= -\eta nc \end{align} Where the physical domain is a simple square, $\Omega = [0,100]\times[0,100]$ and the time domain is $T = [0,t]$. The paper I'm working from says the no-flux boundary condition is: $$\left[D_n\nabla n - \left(\frac{\chi}{1+\alpha c}n\nabla c\right) - (pn\nabla f)\right]\cdot \hat{N} = 0$$ Where $\hat{N}$ is the outward unit normal vector.

What I've done: I've set up the central difference versions of the first 3 equations (Overleaf PDF here), but I'm not sure how to set up $\frac{\partial n}{\partial t}$ on the boundary using the BC above. I've tried converting the BC into a finite difference version and using it to cancel some terms, but I'm not sure if this is the right way to go about this. Any advice is appreciated!


A finite volume viewpoint would be very covenient for your problem. It is not too different from finite difference. Sticking to just 1-d, your first cell would be the interval $[0,h]$ and you store the solution at the center of this cell, say $n_0$. Then the semi-discrete scheme for this cell would be $$ \frac{d n_0}{dt} + \frac{F_{1/2} - F_{-1/2}}{h} + ... = 0 $$ From your boundary condition $F_{-1/2} = 0$. The other fluxes can be computed by a finite difference, e.g., $$ F_{1/2} = D_n \frac{n_1 - n_0}{h} - \frac{\chi n_{1/2}}{1+\alpha c_{1/2}} \frac{c_1 - c_0}{h} - p n_{1/2} \frac{f_1 - f_0}{h} $$ where $n_{1/2}=(n_0 + n_1)/2$, etc.

  • $\begingroup$ $f_{1/2}$ is already given by the ODE though; like Wolfgang mentioned, $f$ and $c$ don't have any spatial derivatives, so this doesn't provide any new equations. I'm trying to estimate $\partial n/\partial t$ which has spatial derivatives and so we need $n_{-1/2}$. I've also thought about doing forward/backward differences on the boundary but wasn't sure if this is recommended. $\endgroup$ – letslearnmath Apr 2 '18 at 18:16
  • $\begingroup$ I should not have used f for flux, I corrected it. You dont need $n_{-1/2}$ because the flux across the boundary face is zero. $\endgroup$ – cfdlab Apr 3 '18 at 2:50

The equations for $f$ and $c$ do not have any spatial derivatives, and consequently they do not require (and cannot have, in fact) any boundary conditions. Thus, only $n$ has boundary conditions, and these are in your case of Neumann type. It is convenient to rewrite the equation you have as $$ \left(\frac{\chi}{1+\alpha c}n\nabla c\right)\cdot\hat{N} = \left[D_n\nabla n - (pn\nabla f)\right]\cdot \hat{N} $$ to make clear what condition has to be satisfied.

In a finite element context, the term on the left hand side of this appears naturally as a result of integration by parts, and dealing with the boundary condition is therefore a trivial matter of substituting the term by what you have here from the boundary condition. I am not enough of a finite difference person to say how best to do this in that context.

  • $\begingroup$ Thank you for your input! I have to use the finite difference method but I've been learning more about the FEM in another class and it's good to think about as well! $\endgroup$ – letslearnmath Apr 2 '18 at 1:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.