# How to apply Neumann boundary condition to wave equation using finite differeces

I am trying to solve 1D wave equation using explicit finite difference scheme and I need to apply Neumann boundary conditions to it:

$$u_{xx} = c^2u_{tt}$$

With initial conditions:

$$u_t(x, t=0) = 0$$

$$u(x, t=0) = 0$$

and boundary conditions:

$$u(x=L, t) = 0$$

$$u_x(x=0, t) = g(t)$$

The finite difference scheme is as follows:

$$u_i ^{(k+1)} = 2u_i ^k - u_i ^{(k-1)} + \alpha ^2(u_{i + 1} ^k - 2u_i ^k + u_{i - 1} ^k)$$

where $$\alpha ^2 = (\Delta t /c \Delta x )^2$$.

With "i" the index for "x" and "k" for "t". Approximating the Neumann BC we have:

$$u_{-1} ^k = u_0 ^k - \Delta x g(t_k)$$

My code for solving this problem is as follows:

def fd_neumann(x, t, c, gt):
dx = x[1] - x[0]
dt = t[1] - t[0]
k2 = (dt/(c*dx))**2.0
assert k2 <= 1.0

Nx, Nt = len(x), len(t)

u = np.zeros((Nx, Nt), dtype = np.float64)

for k in range(1, Nt - 1):
gk = gt(t[k])
u[0, k + 1] = 2.0*u[0, k] - u[0, k - 1] + k2*(u[1, k] - u[0, k] - dx*gk)
for i in range(1, Nx - 2):
u[i, k + 1] = 2.0*u[i, k] - u[i, k - 1] + k2*(u[i + 1, k] - 2.0*u[i, k] + u[i - 1, k])

return u


The result is completly wrong and I just can't find what I am missing.

Normally one uses $$c$$ to represent the wave velocity so you would have $$c^2u_{xx}=u_{tt}$$ but I don't think that's your problem, after all it's just a number. Your big problems seem to be with the boundary conditions. If we start with
$$u_1=u_0+hu_0^{\prime}+\frac12h^2u_0^{\prime\prime}+\frac16h^3u_0^{\prime\prime\prime}+O(h^4)$$
$$u_2=u_0+2hu_0^{\prime}+2h^2u_0^{\prime\prime}+\frac43h^3u_0^{\prime\prime\prime}+O(h^4)$$
We would get $$u_0^{\prime\prime}=\frac{8u_1-u_2-7u_0-6hu_0^{\prime}}{2h^2}+O(h^2)=\frac{8u_1-u_2-7u_0-6hg(t)}{2h^2}+O(h^2)$$ Even a first-order approximation would be $$u_0^{\prime\prime}=\frac{2u_1-2u_0-2hu_0^{\prime}}{h^2}+O(h)=\frac{2u_1-2u_0-2hg(t)}{h^2}+O(h)$$ So you seem to be off by a factor of $$2$$. So the first thing you could try fixing is this boundary condition, either via the first- or second-order method. I would prefer second-order because all your other difference equations are second-order accurate.

I tried coding this up in Matlab:

% wave.m

clear all;
close all;
x0 = 0;
xf = 1;
umin = -1.5;
umax = 1.5;
t0 = 0;
tf = 10;
f = 2;
omega = 2*pi*f;
g = @(t) sin(omega*t);
c = 5;
Nx = 100;
x = linspace(x0,xf,Nx);
dx = (xf-x0)/(Nx-1);
dt = dx/c/1.1246; % stable
%dt = dx/c/1.1245; % unstable
t = t0;
u = zeros(size(x));
plot(x,u);
axis([x0 xf umin umax]);
drawnow;
v = zeros(size(x));
t = t+dt;
v(1) = (4*v(2)-v(3)-2*dx*g(t))/3;
plot(x,v);
axis([x0 xf umin umax]);
drawnow;
while t<tf,
u(end) = 0;
u(1) = 2*v(1)-u(1)+(c*dt/dx)^2*(8*v(2)-v(3)-7*v(1)-6*dx*g(t))/2;
u(2:end-1) = 2*v(2:end-1)-u(2:end-1)+(c*dt/dx)^2*(v(1:end-2)-2*v(2:end-1)+v(3:end));
plot(x,u);
axis([x0 xf umin umax]);
t = t+dt;
drawnow;
v(end) = 0;
v(1) = 2*u(1)-v(1)+(c*dt/dx)^2*(8*u(2)-u(3)-7*u(1)-6*dx*g(t))/2;
v(2:end-1) = 2*u(2:end-1)-v(2:end-1)+(c*dt/dx)^2*(u(1:end-2)-2*u(2:end-1)+u(3:end));
plot(x,v);
axis([x0 xf umin umax]);
t = t+dt;
drawnow;
end


And it seemed to become von Neumann unstable when the equivalent of your k2 variable was $$1/1.1245^2$$ I'm not sure just why it happened at that point because it should be stable at all points except the left boundary when $$k2\le1$$ but the left boundary becomes unstable much earlier than this. The solution goes crazy as well if driven in resonance.

Just what is going wrong in your solution? have you tried just decreasing $$\Delta t$$ to see if it's von Neumann unstable, or changing the driving frequency to see if you are in resonance?

• Thank you for your answer, but do you mind commenting your code? I ran it on Matlab and worked pretty fine, but I didn't figure out whats the v term means and why did you use that expression.
– Gabs
Commented Mar 10, 2020 at 19:51
• My code is pretty much the same as yours except for the handling of the Neumann boundary condition. Rather than keeping the data for all previous times around when we're just gonna do an animation, I just use the $2$ most recent steps because they are all that are needed for the next step. The first step is $u$, the second is $v$, the third is $u$ again and so on. I overwrite the data of the first step when computing the third step, but I don't need it any more. Did you try uncommenting the line marked %unstable? It's really cool how going just over the edge like this makes every thing go bl Commented Mar 10, 2020 at 20:26
• Also did you make some progress in diagnosing your python code? Have you tried changing the boundary difference equation, time step, and driving frequency? It's kind of hard on this end to see what is going wrong without being able to see the whole thing. Commented Mar 10, 2020 at 20:30
• OK, but can you describe in details how did you treat the Neumann boundary conditions?
– Gabs
Commented Mar 12, 2020 at 0:45
• Well. I thought I did. I wrote down the Taylor polynomial of degree $3$ for $u_1=u(x_1,t_n)$ and $u_2=u(x_2,t_n)$ and then solved for $$u_0^{\prime\prime}=\left.\frac{\partial^2u}{\partial x^2}\right|_{x=x_0,t=t_n}$$ in terms of $u_0$, $u_1$, $u_2$, and $u_0^{\prime}$, eliminating $u_0^{\prime\prime\prime}$. Then I could substitute the Neumann boundary condition $u_0^{\prime}=\left.\frac{\partial u}{\partial x}\right|_{x=x_0.t=t_n}=g(t_n)$. Then I had an expression for $u_0^{\prime\prime}$ that I used just like the expression for $u_i^{\prime\prime}$ Commented Mar 12, 2020 at 0:59