# Is this a valid way to implement Neumann BCs in finite differences?

I'm trying to solve the 1D heat equation with Neumann boundry conditions numerically using finite differences: $$u_t = \alpha u_{xx}$$ $$u_{x}(0, t) = u_{x}(L, t) = k$$ $$u(x, 0) = u_0(x)$$ The main problem here is that $$u_x$$ isn't in the equation, so I thought I could use the Taylor expansions of $$u(x)$$: $$u(x+h) = u(x) + u_x(x)h + u_{xx}(x)\frac{h^2}{2} + u_{xxx}(x)\frac{h^3}{6}$$ $$u(x+2h) = u(x) + 2u_x(x)h + 2u_{xx}(x)h^2 + u_{xxx}(x)\frac{4h^3}{3}$$ And then I cancel out the $$h^3$$ terms by doing: $$u(x + 2h) - 8u(x + h) = -7u(x) - 6u_x(x)h - 2u_{xx}(x)h^2$$ Then solving for $$u_{xx}$$ yields: $$u_{xx}(x) = \frac{-7u(x) + 8u(x+h) - u(x+2h) - 6u_x(x)h}{2h^2}$$ And since we know $$u_x(0) = k$$ (at the boundry), we get the formula: $$u_{xx}(0) = \frac{-7u(0) + 8u(h) - u(2h) - 6kh}{2h^2}$$ A similar formula can be found for $$u_x(L)$$. The reason I wanted to cancel out the $$h^3$$ terms is because I wanted a formula with and error proportional to $$O(h^4) / h^2 = O(h^2)$$ since it's what's used on the interior points.

So now to the question: is this a valid way of incorporating the Neumann BC or is it unstable in any serious way?

• The way you derived that approximation of the second derivative is quite unusual for me, I've never seen it before. Euristically, the lack of simmetry in your approximation could be a problem, but don't trust me. Anyway, I posted an answer explaining how I'd do – VoB Jan 13 '20 at 18:49
• I have not seen this type of scheme. But all your steps are correct. Write the semi-discrete FD scheme as $u_t = A u + b$ where $b$ contains the boundary condition terms. If the eigenvalues of $A$ are real and negative, then you have stability. To begin, you can numerically check the eigenvalues on some grid. And if it indicates stability, perhaps you can try to prove this analytically. – cfdlab Jan 16 '20 at 6:11
• Have you tried any numerical testing? – Gabriele Jan 22 '20 at 16:50

For any time you want $$\frac{u(t,L+h) - u(t,L-h)}{2h} = K , \quad \star$$Of course $$x=L+h$$ is not a point in your computational domain, but from the equation above you can get an equation for $$u(x, L+ h)$$ and substitute in.

More precisely, if you discretize in space (with finite difference) from $$x_0, \ldots, x_n$$, then $$x_{n+1} = L+h$$ and $$x_{-1} = L-h$$.

Now, using a second order approximation of space derivative you have $$u_{xx}(t,x_i) = \frac{u_{i+1} - 2 u_i + u_{i-1}}{h^2}$$.

Of course problems arise when you're at the boundaries. Let's see what we can do:

When you are at $$x=x_0$$ you have:

$$u_{xx}(x_0) = \frac{u_1-2u_0 + u_{-1}}{h^2}, \quad i =0, \ldots, N$$

But $$u_{-1}$$ can be computed from $$\star$$ using the information on the boundary derivative.

The very same argument has to be done for $$x_{N}$$

Now you can set up your matrix, where only the first and last row changes and integrate in time with a suitable numerical method (notice that semidiscretized problems are usually stiff, hence explicit integrators could be quite expensive)

• Thank you for your answer, I will test this method out as well :-) – Hugo Gransträm Jan 14 '20 at 18:37
• But I don't know about you, but I don't think this answers the question I posed. But it does offer an alternative way that probably is better :-) Thank you for taking your time ;-) – Hugo Gransträm Jan 14 '20 at 18:39
• My method is the one you would see in a numerical analysis class. A possibility to implement Neumann bdary conditions could also just use a one-sided approximation of the derivative at the boundaries, namely forcing $u_0 - u_1 =0$ and $u_N - u_{N-1} = 0$, which is a first order approximation – VoB Jan 17 '20 at 18:04
• So your's definitely a better method ;-) – Hugo Gransträm Jan 18 '20 at 17:31