# Neumann boundary conditions diffusion equations methods of lines

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

• 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. – boyfarrell Nov 20 '17 at 21:57

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.

• Trapezoidal rule isn't required for this to work. – Chris Rackauckas Nov 21 '17 at 4:00
• 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 – alessio lapolla Nov 21 '17 at 16:52
• 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. – cfdlab Nov 22 '17 at 2:33