Numerical method with a time derivative boundary condition

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.

Cheers!

• 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? – Phill May 26 '15 at 2:30
• Also, can you explain why this is a Dirichlet condition at $x=0$? – Phill May 26 '15 at 2:36
• Strictly speaking, it's not, but since it olnly involves the value of the function at $x=0$, I decided to call it so! – Dmoreno May 26 '15 at 2:43
• Fair enough, I guessed it was something like that :D – Phill May 26 '15 at 2:56
• Happy to help @Phill! – Dmoreno May 26 '15 at 3:28