Neumann boundary conditions diffusion equations methods of lines

Question

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

Answer 1

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.

Answer 2

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.