Using Finite Difference method for 1d diffusion equation

Question

In the book "A Primer on PDEs" by Salsa et al. they use the finite difference method to solve the one dimensional non-homogeneous diffusion equation $$ u_t(x,t)-u_{xx}(x,t)=f(x,t) $$ For $x \in \, (0,1)$ with Dirichlet boundary condition $u(0,t)=0$ and Neumann boundary condition $u_x(1,t)=0$ and a smooth initial condition. Then the book introduces an $N$ component vector $\mathbf{U}(t)$ such that $U_i(t)=u(x_i,t)$ where $x_i$ is the discretized spatial coordinate and $1\leq i\leq N$. For the Neumann boundary condition the book uses the following approximation for the first order derivative at the boundary point ($i=N$): $$u_x(x_N,t)=\frac{1}{2h}(-3u(x_{N-2},t)+2u(x_{N-1},t)-u(x_N,t))$$ And writes the discretized equation as: \begin{equation} \dot{\mathbf{U}}(t)+\mathbf{A}\mathbf{U}(t)=\mathbf{F}(t) \,\,\, \, \, \, \, \, \, \, \, \, (i) \end{equation} Where $$\mathbf{A}=\frac{1}{\Delta x^2}\begin{bmatrix} 2 & -1 & 0 & \dots & \dots & \dots & 0 \\ -1 & 2 & -1 & 0 & \dots & \dots & 0\\ 0 & -1 & 2 & -1 & 0 & \dots & 0\\ \vdots & & \ddots & \ddots & \ddots & & \vdots \\ 0 & \dots & & -1 & 2 & -1 & 0\\ 0 & \dots & & 0 & -1 & 2 & -1 \\ 0 & \dots & & 0 & -\frac32 & 2 & -\frac12 \end{bmatrix}$$ Where $\Delta x = \frac{1}{N}$ and $$ \mathbf{F}(t)=\begin{pmatrix} f(x_1,t) \\ f(x_2,t) \\ \vdots \\ f(x_{N-1},t) \\ 0 \end{pmatrix}$$

My problem is with that last row of $\mathbf{A}$. I know that it has something to with the Neumann boundary condition, but it doesn't imply that the Neumann condition gets satisfied. If we put $\mathbf{A}$ into $(i)$ we find a relation for the time derivative of $u(x_N,t)$. So I would like to know that how the Neumann condition is satistied in equation $(i)$.

Answer 1

The last row of the matrix gives $$-\frac{3}{2}u_t(x_{N-2})+2u_t(x_{N-1})-\frac{1}{2}u_t(x_N)$$ which is a second-order expression for $(u_x)_t(x_N)$, and the last entry of $\mathbf{F}$ is 0. This is a discrete statement that $(u_x)_t=0$ at $x_N$. You are correct that there is something fishy about this. It does not work unless $u_x=0$ initially, which perhaps it is.

A different way to do it would be to ignore the last row and compute $u_t$ and hence $u$ for $n=0,1..N-2,N-1$, and then find $u_N$ by extrapolation

It is very hard to write a textbook with no errors at all.

Answer 2

The way I was taught to think about a Neumann condition in the finite difference framework is to insert a "ghost point" one mesh point outside the grid in the normal direction. This ghost point then contributes to the system in two places. First, it contributes to the Laplacian at the point on the "real boundary" of the region ($x_N$ in this case). Second, one can (and should) choose the discretization of the Neumann condition so that the ghost point contributes to the Neumann condition; this determines the value of the solution at the ghost point (although as the name implies, we do not really care what this value is).

In this example problem, if we take the centered difference for the Neumann condition, then it is requiring $\frac{u_{N+1}-u_{N-1}}{2h}=0$ so $u_{N+1}=u_{N-1}$. One way to set up the method is to not solve for $u_{N+1}$ and instead replace any appearances of it in the system with $u_{N-1}$. In this problem that means that the Laplacian equation $\frac{du_N}{dt}-(u_{N-1}-2u_N+u_{N+1})=h^2 f(x_N,t)$ becomes instead $\frac{du_N}{dt}-(2u_{N-1}-2u_N)=h^2 f(x_N,t)$.

This ghost point concept is closer to how finite element/finite volume methods work, and does not require anything of the initial data. By contrast, if we do not "force" things like this then the given initial data may violate the Neumann condition, and then problems can arise as you seem to have noticed. In this situation our job as numerical programmers becomes a headache, because we have to inspect the initial data to determine whether it is suitable and yell at the user if it is not.