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.

  • $\begingroup$ Equation [7] looks fine to me. $\endgroup$ Dec 29, 2015 at 22:20
  • $\begingroup$ Then you would first solve eq. [7] and then eq. [6] to obtain the derivative and N-1? $\endgroup$
    – Toulousain
    Dec 29, 2015 at 23:21
  • $\begingroup$ [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. $\endgroup$ Dec 30, 2015 at 2:37
  • $\begingroup$ 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 $\endgroup$
    – Iftekhar
    Aug 3, 2017 at 8:24

1 Answer 1


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] $$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.