# 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[1]$
$c(x,0)=1 \quad[2], \quad \text{ initial condition}$
${dc\over dx} \mid_{(0,t)} =0\ \text{ for t>0,} \:[3], \text{ symmetry at the center}$
$D{dc\over dx} \mid_{(L,t)} =k_L(C_L-C\mid_{(L,t)}) \quad[4]$

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[5]$
${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[6]$

And then, for the last condition I am not sure how to implement it. Discretizing equation [4], 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[7]$

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 [7] looks fine to me. Dec 29, 2015 at 22:20
• Then you would first solve eq. [7] and then eq. [6] to obtain the derivative and N-1? Dec 29, 2015 at 23:21
• [6] is continuous in time, so you still have some discretizing to do. You can view the full discretization as a coupled system or use [7] to eliminate an unknown. Dec 30, 2015 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 Aug 3, 2017 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 [6] is correct, I suggest to change [5] and [7] 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 [5]$$ and $$\frac{dc_N}{dt}=D \frac{-2c_N+2c_{N-1}}{dx^2}+\frac{2k_L}{dx}(C_L-c_N)+r \quad [7]$$