I want to solve the diffusion equation using the method of lines with Neumann boundary conditions $$ \frac{\partial p}{\partial t}=\frac{\partial^2p}{\partial x^2}\\ \frac{\partial p}{\partial x}(x=0)=\frac{\partial p}{\partial x}(x=1)=0\\ p(x,0)=\begin{cases} 1/dx \quad x=0.5,\\ 0 \quad elsewhere \end{cases} $$ For sake of simplicity I use the forward Euler so I choose $dx=0.1, \; dt=0.001$ to ensure the stability. Using a second order scheme I solve the boundary conditions imposing $$ \frac{p_1-p_{-1}}{2 dx}=0,\\ \frac{p_{N+1}-p_{N-1}}{2 dx}=0.\\ $$ So the tridiagonal matrix looks like $$ \frac{1}{dx^2} \begin{pmatrix} -2&2&0&0\\ 1&-2&1&0\\ 0&1&-2&1\\ 0&0&2&-2 \end{pmatrix}. $$ At this point the solution of my equation is given by $$ \vec{p}^{k+1}=dt(T\vec{p}^k)+\vec{p}^k. $$ This boundary conditions impose the conservation of probability but it seems that is not the case since for different $N$ (number of grid points on the $x$ axis) I obtain different probabilities. For example, for $N=10,20,100,1000$ I obtained that the total probability is $1.11111,1.05263,1.00012,1.0000$. It seems that I need a really fine grid to obtain the correct physical result.

Did I make a conceptual mistake or I just made an error in my code? Thank for the answers

  • $\begingroup$ Your discretised system is an approximation and the finite difference method is not exactly conservative. If use used finite volume method, which is conservative discretisation it would improve. I think 1000 points is not unreasonable. $\endgroup$
    – boyfarrell
    Nov 20, 2017 at 21:57

2 Answers 2


It is easy to see that the function \begin{equation} P(t) := \int \limits_0^1 p(x,t) \, \mathrm{d}x, \quad t \geq 0, \end{equation} is constant over time (i. e. $P(t) = P(0)$, $\forall \, t \geq 0$), if $p$ is a solution of your initial-boundary-value problem.

If you discretize the problem in space as you described, then you obtain a system of linear ordinary differential equations in time of the form \begin{equation} \frac{d\boldsymbol{p}}{dt} = \boldsymbol{\underline{T}} \boldsymbol{p} \end{equation} for the vector-valued function $\boldsymbol{p}(t) \in \mathbb{R}^{N+1}$, with the tridiagonal matrix $\boldsymbol{\underline{T}} \in \mathbb{R}^{(N+1) \times (N+1)}$ that you have specified.

Now if you define some function $\tilde{P}(t) := \boldsymbol{c}^{\top} \boldsymbol{p}(t)$, $t \geq 0$, with a vector $\boldsymbol{c} \in \mathbb{R}^{N+1}$, then this function will be constant over time (i. e. $\tilde{P}(t) = \tilde{P}(0)$, $\forall \, t \geq 0$) if and only if $\frac{d\tilde{P}}{dt} = 0$. The time derivative of $\tilde{P}$ is given by \begin{equation} \frac{d\tilde{P}}{dt} = \frac{d}{dt}(\boldsymbol{c}^{\top}\boldsymbol{p}) = \boldsymbol{c}^{\top} \frac{d\boldsymbol{p}}{dt} = \boldsymbol{c}^{\top} \boldsymbol{\underline{T}} \boldsymbol{p}, \end{equation} and the right-hand side is zero for any $\boldsymbol{p}$ if and only if $\boldsymbol{c}^{\top} \boldsymbol{\underline{T}} = \boldsymbol{0}_{\mathbb{R}^{N+1}}^{\top}$. For your tridiagonal matrix $\boldsymbol{\underline{T}}$ this is satisfied if and only if \begin{equation} \boldsymbol{c}^{\top} = \left( \begin{array}{ccccc} \frac{1}{2} & 1 & \cdots & 1 & \frac{1}{2} \end{array} \right) \quad \Rightarrow \quad \tilde{P} = \boldsymbol{c}^{\top} \boldsymbol{p} = \frac{1}{2} p_0 + \sum_{i=1}^{N-1} p_i + \frac{1}{2} p_N. \end{equation} This corresponds to using the trapezoidal rule for the approximation of the integral in $P$. If you calculate your $\tilde{P}$ in any other way then its value will not be constant over time for this discretization.


Consider a uniform grid $$ x_i = ih, \quad i=0,1,\ldots,N $$ Your equations are $$ dp_0/dt = 2(p_1 - p_0)/h^2 $$ $$ dp_i/dt = (p_{i-1} - 2 p_i + p_{i+1})/h^2, \quad 1 \le i \le N-1 $$ $$ dp_N/dt = 2(p_{N-1} - p_N)/h^2 $$ This can also be written as $$ dp/dt = Tp, \qquad p = [p_0, p_1, \ldots, p_N]^\top $$ where $T$ is tridiagonal matrix. Compute integral of $p$ by Trapezoidal rule $$ I(t) = \left(p_0/2 + \sum_{i=1}^{N-1} p_i + p_N/2\right) h $$ Using the scheme $$ dI/dt = \frac{1}{h} \left(p_1 - p_0 + \sum_{i=1}^{N-1} (p_{i-1} - 2 p_i + p_{i+1}) + p_{N-1} - p_N \right) $$ Rearrange the sum to show that $$ dI/dt = 0 $$

It is not clear from the question how the total probability was calculated. As shown above, using Trapezoidal rule, the numerical scheme will conserve the total probability.

  • $\begingroup$ Trapezoidal rule isn't required for this to work. $\endgroup$ Nov 21, 2017 at 4:00
  • $\begingroup$ To calculate the the probability I just summed the entries of the vector p after a fixed time of time steps and then I multiplied by dx to take in account the normalization. I totally do not understand what you integrate; I do not understand how you find the solution of the PDE $\endgroup$
    – apelle
    Nov 21, 2017 at 16:52
  • 1
    $\begingroup$ I have written each equation from your matrix form. You are calculating total probability wrongly. Do it with Trapezoidal rule and your scheme will conserve total probability. $\endgroup$
    – cfdlab
    Nov 22, 2017 at 2:33

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.