# Using Finite Difference method for 1d diffusion equation

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: $$\dot{\mathbf{U}}(t)+\mathbf{A}\mathbf{U}(t)=\mathbf{F}(t) \,\,\, \, \, \, \, \, \, \, \, \, (i)$$ 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)$.

• When $i=N$, this means $x_N=1$ and the assumption is $u_x(1,t)=0$, thus $u_x(x_N,t)=0$. Now you can look at the last equation in the system, which corresponds to the last row of the matrix. – Frank Lu Jul 5 '17 at 3:54
• The last equation gives $u_t(1,t)$ and doesn't imply that $u_x(1,t)=0$. – Hosein Jul 5 '17 at 4:02

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)$.
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