# Using Neumann boundary conditions for the wave equation

I have the following code to solve the wave equation in 2D:

range=2;
trange=5;
region=Region[Rectangle[range{-1,-1},range{1,1}]];
wave2D = NDSolveValue[
{D[u[t,x,y] ,{t,2}]-
Laplacian[u[t, x, y], {x, y}] == 0,
u[0, x, y] == E^(-5 (x^2 + y^2)),
Derivative[1, 0, 0][u][0, x, y] == 0
(* absorbing condition *)
},
u, {t, 0, trange}, {x, y} \[Element] region] // Quiet;


I want to add an absorbing boundary condition: $$\partial_t u=\partial_x u$$ on the boundary and I think I have to use NeumannValue somehow but I don't understand the documentation of NeumannValue. Can someone show me how to do that and explain the idea behind NeumannValue in simple terms?

Edit: with the help of @Ulrich Neumann (perfect name for this question) I was able to make the following plot

plot = Table[
Plot3D[wave2D[t, x, y], {x, -range, range}, {y, -range, range},
PlotRange -> {-.5, 1}, ImageSize -> 300], {t, 0, 3.8, .075}];
ListAnimate[plot]


• The acoustics tutorials have ton of information on this. Have look there. Apr 30 at 10:52
• What should "Derivative[1, 0, 0][u] == 0" mean? A function needs arguments. Apr 30 at 11:07
• Have a look here or here Apr 30 at 11:11
• @DanielHuber Oops, I edited it now. Apr 30 at 14:24
• Related: mathematica.stackexchange.com/a/245309/1871 As to ABC, see also: mathematica.stackexchange.com/a/236112/1871 May 1 at 4:07

With NeumannValue the absorbing boundary conditions might be formulated for x==2and x==-2. Roughly NeumannValue describes the flux perpendicular to the boundary:

    range = 2;
trange = 5;
region = Region[Rectangle[range {-1, -1}, range {1, 1}]];
wave2D = NDSolveValue[{D[u[t, x, y], {t, 2}] -
Laplacian[u[t, x, y], {x, y}] ==
NeumannValue[-Derivative[1, 0, 0][u][t, x, y], x == 2]+NeumannValue[ -Derivative[1, 0, 0][u][t, x, y], x == -2],
u[0, x, y] == E^(-5 (x^2 + y^2)),
Derivative[1, 0, 0][u][0, x, y] == 0
(*absorbing condition*)},
u, {t, 0, trange}, {x, y} \[Element] region]

Manipulate[Plot3D[wave2D[t, x, y], Element[{x, y}, region],PlotRange -> {-1, 1}], {{t, 2 }, 0, 5, Appearance -> "Labeled"}]


• Very nice! I want it along all boundaries but that works by setting the predicate to true. What I still think is weird is that you don't need the term laplacian u == 0 anymore, does NeumannValue evaluate to zero for the region that is not on the boundary? Apr 30 at 14:29
• @AccidentalTaylorExpansion Laplacian[u[t, x, y], {x, y}]is still part of the pde! Apr 30 at 14:31
• I meant to say $\partial_{tt}u-\nabla^2 u=0$ does not occur on its own anymore, only in conjunction with NeumannValue, which would imply for me that NeumannValue evaluates to zero inside the bulk and to something else (in this case $\partial_t u$ on the boundary). It's still a bit confusing how it works but I guess it works now Apr 30 at 15:06

As mentioned in the comment above, the problem is related with (and essentially equivalent to, in my view) this one, but let me elaborate with an answer anyway.

The explanation for the usage of NeumannValue can be found in Details section of document of NeumannValue, and the tutorial Finite Element Method Usage Tips (Particularly Formal Partial Differential Equations and NeumannValue and Formal Partial Differential Equations section) but maybeprobably a bit hard to follow. In short, the NeumannValue is defined based on a so-called formal PDE, which can be checked with the new-in-12.2 NDSolveFEMGetInactivePDE. (If you're not yet in v12.2, try the function in this post. )

Let's check the formal PDE corresponding to the 2D wave equation:

<< NDSolveFEM;

{eqLHS, ic} =
With[{u = u[t, x, y]}, {D[u, {t, 2}] -
Laplacian[u, {x, y}], {u == E^(-5 (x^2 + y^2)), D[u, t] == 0} /. t -> 0}];

GetInactivePDE@
First@NDSolveProcessEquations[{eqLHS == 0, ic},
u, {t, 0, trange}, {x, y} ∈ region]


As mentioned in Details section of document of NeumannValue:

$$∇·(-c ∇u-α u+γ)+…=f+\text{NeumannValue}[g-q u,\text{pred}]$$ is used to specify the flux over the part of the boundary $$∂Ω$$ where pred is true, such that $$\overset{\rightharpoonup }{n}·(c ∇u+α u-γ)=g-q u$$ holds.

Clearly, in our case $$-c ∇u-α u+γ=$$-{{1, 0}, {0, 1}}.Inactive[Grad][u[t, x, y], {x, y}].

Next step is to determine the $$g-q u$$. For $$x=2$$, the needed 1st order ABC is $$\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=0$$, $$\overset{\rightharpoonup }{n}=(1,0)$$, so

{1, 0} . (-(-{{1, 0}, {0, 1}} . Inactive[Grad][u[t, x, y], {x, y}])) == gMinuspu


should be equivalent to

With[{u = u[t, x, y]}, D[u, t] + D[u, x] == 0]


at $$x=2$$. Then it's clear gminuspu i.e. $$g-pu$$ should be -D[u[t, x, y], t] at $$x=2$$.

Similarly we'll find at all the boundary of the rectangle region, the needed $$g-pu$$ terms are all -D[u[t, x, y], t]. So the needed NeumannValue is

sol = NDSolveValue[{eqLHS == NeumannValue[-D[u[t, x, y], t], True], ic},
u, {t, 0, trange}, {x, y} ∈ region];

plot = Table[
DensityPlot[sol[t, x, y], {x, -range, range}, {y, -range, range},
PlotRange -> {-.5, 1}, ImageSize -> 300, PlotPoints -> 50], {t, 0, 3.8, .075}];
ListAnimate[plot]


From the animation we can see the reflection is still somewhat obvious, this is expected because we're using 1st order ABC. To reduce the reflection, you may consider 2nd order ABC as shown here. Setting up PML as shown in the tutorial Acoustics in the Time Domain is also a possible choice.

• (+1) do you have suggestions for how to improve the documentation to get rid of the 'probaby hard to follow'? May 3 at 16:53
• @user21 One thing I can think of is, perhaps the document should emphasize a bit more about the importance of formal PDE/Inactive PDE. Adding example showing how to determine the NeumannValue with the help of GetInactivePDE (for certain strange PDE) may improve the readability of document, too. As we've seen, checking the corresponding formal PDE is almost unavoidable for setting NeumannValue` correctly. May 3 at 17:04
• Thanks, I added some more information to the documentation. I feel that at the end of the day one has to think through the derivation of the FEM like it is hinted add in the NeumannValue ref page or the FEM tutorials to really understand the meaning of NeumannValue and how to use that. I am not sure there is really a short cut. But again, I tried to add more information to the documentation and made GetInactivePDE more prominent in conjunction with the sections on understanding PDE parsing and NeumannValue. May 4 at 5:32