# Implementation of Neumann boundary condition with method of lines - 1D diffusion/reaction equation

I am solving the monodimensional diffusion/reaction equation by discretization using the well-known method of lines

${\partial c\over\partial t}=D{\partial^2 c\over\partial x^2}+ r\text{ for t>0 and -L<x<L} \quad$
$c(x,0)=1 \quad, \quad \text{ initial condition}$
${dc\over dx} \mid_{(0,t)} =0\ \text{ for t>0,} \:, \text{ symmetry at the center}$
$D{dc\over dx} \mid_{(L,t)} =k_L(C_L-C\mid_{(L,t)}) \quad$

This last equation represents that the diffusion at the boundary equals the external mass transfer but I am having trouble to transform it to the method of lines. The discretization of the system becomes:
${dc_i\over dt} =D{2 c_{i+1}-c_i\over \Delta x^2}+r \text{ for i=1} \quad$
${dc_i\over dt} =D{ c_{i+1}-2c_i+c_{i-1}\over \Delta x^2} +r\text{ for i = 2...N-1} \quad$

And then, for the last condition I am not sure how to implement it. Discretizing equation , I obtain $C_i ={k_L\:c_L + {D\over \Delta x} c_{i-1} \over k_L + {D\over \Delta x} } \text{ for i = N} \quad$

I don't think that this is correct given that $C_N$ does not depend on time, is calculated indirectly at each time step from $C_{N-1}$. And that I cannot obtain the derivative for ${dc_{N-1}\over dt}$. But I cannot see how that condition should be discretized. Thank you for your help.

• Equation  looks fine to me. – David Ketcheson Dec 29 '15 at 22:20
• Then you would first solve eq.  and then eq.  to obtain the derivative and N-1? – Toulousain Dec 29 '15 at 23:21
•  is continuous in time, so you still have some discretizing to do. You can view the full discretization as a coupled system or use  to eliminate an unknown. – David Ketcheson Dec 30 '15 at 2:37
• Thank you very much, Peter, for your answer. Your answer helped me in brushing up my knowledge. I would like to mention one point. When I applied the 2nd Boundary Condition, I acheived this: dc(N)/dt = D[-c(N)+c(N-1)]/dx^2 + kL/dx*[CL-c(N)] + r – Iftekhar Aug 3 '17 at 8:24

## 1 Answer

First a general comment. In the case of Neumann boundary condition, opposite to Dirichlet boundary condition where you specify the value of unknown function, you should approximate the PDE also in boundary nodes, it means you should not try to eliminate the unknown functions there as you suggest.

Furthermore, you should use the central difference to approximate the derivative in Neumann boundary condition, otherwise you lose accuracy that is not necessary.

Although in principle your method can work, you might need a quite fine mesh to obtain a good approximation, because it is like you solve it by MOL on a domain $(h,L-h)$ and with only a first order accurate scheme.

The following scheme is second order accurate. The scheme  is correct, I suggest to change  and  as follows.

Let the index $i$ denote in which space coordinate $x_i$ you approximate your equation. I prefer $x_i = i h$ for $h=L/N$ and $i=0,1,\ldots,N$. Your aim should be to construct a system of ODEs for unknown functions $c_i(t) \approx c(x_i,t)$ for $i=0,1,\ldots,N$, it means including the concentrations at boundary nodes.

Using central difference to approximate $dc/dx$ in Neumann boundary conditions and using it to eliminate auxiliary functions $c_{-1}(t)$ and $c_{N+1}(t)$ you should obtain $$\frac{dc_0}{dt}=D \frac{-2c_0+2c_1}{dx^2}+r \quad $$ and $$\frac{dc_N}{dt}=D \frac{-2c_N+2c_{N-1}}{dx^2}+\frac{2k_L}{dx}(C_L-c_N)+r \quad $$