# Implementing no-flux boundary condition reaction-diffusion PDE

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.

• $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. – letslearnmath Apr 2 '18 at 18:16
• 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. – 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.

• 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! – letslearnmath Apr 2 '18 at 1:38