I'm trying to reproduce a result from a paper I'm reading using a numerical scheme that I'm coding myself. The equation is a reaction diffusion PDE. $$\frac{\partial M}{\partial t}=\frac{\partial^2 M}{\partial x^2}+f(M,R)$$ $$\frac{\partial R}{\partial t}=-f(M,R)$$ Where $M$ and $R$ are concentrations of the diffusing and reacting species respectively, I decided on the Crank-Nicolson method for the diffusion part, and forward euler for the reaction part.

In my test cases with the CN scheme I had Neumann BCs, which I discretised as follows: $$\frac{M^{k}_{-1}-M^{k}_{1}}{2\Delta x}=\alpha \to M^{k}_{-1}=2\Delta x \alpha+M^{k}_{1}$$ (This also got rid of the spatial index that fell outside of the domain)

But now I'm faced with a boundary condition like this: $$\frac{\partial M}{\partial t}(0,t)=c+f(M,R)$$ This is essentially the equation for $M$ with the diffusion term replaced by a constant $c$.

How do I deal with this BC in terms of discretisation and removing the negative index?


Simply use your discretization again to have:

$$ \frac{M^{n+1}_0 - M^n_0}{\Delta t} = c + f^n_0,$$ from wich you can determine $M^{n+1}_0$ as a function of what happened before in the boundary, substituting this value in your finite difference scheme for $i = 0$. For instance, you will have:

$$\frac{M^{n+1}_i - M^n_i}{\Delta t} = \frac{1}{\Delta x^2} (M^{n+1}_{i+1} - 2 M_i^{n+1}+ M_{i-1}^{n+1} + M^{n}_{i+1} - 2 M_i^{n}+ M_{i-1}^n ) + f^n_i ,$$ this equation is valid for $i \geq 1$, since we have a Dirichlet condition at $x = 0$. Notice that when you plug $i = 1$ in the equation above, a $M_0^{n+1}$ term appears, which should be written in terms of $M_0^n$ according to the boundary condition. That gives you an equation for $M_0$ at each $n$.

Writing this in matrix form helps a lot.

Hope this helps.


  • $\begingroup$ Thanks for that! So effectively, instead of doing an $i=0$ case and dealing with the negative index by calling on the Neumann boundary conditions, I only go to $i=1$ and use the time dependant boundary conditions to resolve the $i=0$ case? $\endgroup$ – Phill May 26 '15 at 2:30
  • $\begingroup$ Also, can you explain why this is a Dirichlet condition at $x=0$? $\endgroup$ – Phill May 26 '15 at 2:36
  • $\begingroup$ Strictly speaking, it's not, but since it olnly involves the value of the function at $x=0$, I decided to call it so! $\endgroup$ – Dmoreno May 26 '15 at 2:43
  • 1
    $\begingroup$ Fair enough, I guessed it was something like that :D $\endgroup$ – Phill May 26 '15 at 2:56
  • $\begingroup$ Happy to help @Phill! $\endgroup$ – Dmoreno May 26 '15 at 3:28

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.