Absorbing Boundary Condition for Complex/Coupled PDE

Context: This question is relevant to the physical problem of an excited leaky cavity. The boundary condition on the inside of the cavity is Dirichlet. The other end of the cavity is partially blocked by a potential, and beyond that waves are free to escape to infinity. The dynamics are governed by a modified wave equation.

First set the geometry of the problem, with an example Gaussian perturbation and potential:

{xin, xmax, tmax} = {-300, 400, 3000};
{x0, xp, w, wp, v} = {-150, 0, 2, 8, 1};
c = 0.4;
g[x_,t_] := Exp[-(((x - x0) - v*t)/w)^2]
gp[x_] := Exp[-((x - xp)/wp)^2]

And define the modified wave equation:

weqn = D[psi[x, t], {x, 2}] - D[psi[x, t], {t, 2}] -
2 c*I*D[psi[x, t], {t, 1}] + (c^2 - gp[x])*
psi[x, t] == -NeumannValue[-Derivative[0, 1][psi][x, t],
x == xmax];

Because waves should be purely outgoing outside of the potential, to avoid integrating out to large x we can implement an absorbing boundary condition by putting in the NeumannValue[...] term above for some xmax sufficiently far from the potential.

Initial conditions:

ic = {psi[x, 0] == g[x, 0],
Derivative[0, 1][psi][x, 0] == Derivative[0, 1][g][x, 0]};

Boundary condition for the inside of the cavity:

bcd = {psi[xin, t] == 0};

Attempt at a solution:

The absorbing boundary condition seems to require a FEM method, so we write

sol = psi /.
NDSolve[{weqn, ic, bcd}, psi, {x, xin, xmax}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {psi -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 1}}}][[1, 1]];

This works when c = 0, for which weqn is a simple wave equation with a potential. But when c = 0.4 the following warnings are generated,

NDSolve::femdpop: The NDSolveFEMFEMStiffnessElements operator failed.

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

The problem seems to be due to the third (imaginary) term in the PDE. So a natural fix would be to split the equation into its real and imaginary parts and solve the coupled PDE, so that all terms are real. In other words, we take psi = psiR + I*psiI and thus obtain

weqn1 = D[psiR[x, t], {x, 2}] - D[psiR[x, t], {t, 2}] +
2*c*D[psiI[x, t], {t, 1}] + (c^2 - gp[x])*
psiR[x, t] == -NeumannValue[-Derivative[0, 1][psiR][x, t],
x == xmax];
weqn2 = D[psiI[x, t], {x, 2}] - D[psiI[x, t], {t, 2}] -
2*c*D[psiR[x, t], {t, 1}] + (c^2 - gp[x])*
psiI[x, t] == -NeumannValue[-Derivative[0, 1][psiI][x, t],
x == xmax];
ic1 = {psiR[x, 0] == g[x, 0],
Derivative[0, 1][psiR][x, 0] == Derivative[0, 1][g][x, 0]};
ic2 = {psiI[x, 0] == 0, Derivative[0, 1][psiI][x, 0] == 0};
bcd1 = {psiR[xin, t] == 0};
bcd2 = {psiI[xin, t] == 0};

Now we solve

solP = NDSolve[{weqn1, weqn2, ic1, ic2, bcd1, bcd2}, {psiR, psiI}, {x,
xin, xmax}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {psiR -> 2, psiI -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 1}}}][];

Which gives no errors. Then we combine to get the full solution

solR = psiR /. solP[];
solI = psiI /. solP[];
sol[x_,t_] := solR[x,t] + I*solI[x,t];

The problem is, now the "absorbing" boundary condition at xmax is partially reflecting, as can be seen in the animation

Animate[Plot[Abs@sol[x,t], {x, xin, xmax},
PlotRange -> {All, {0, 0.1}}], {t, 0, 800}]

The question is, how should an absorbing boundary condition be properly implemented for this case of coupled PDEs?

Edit: To clarify, the problem here is for a leaky cavity with the above PDE on the inside of the potential, but with the simple wave equation on the outside.

• Why do you define absorb b.c. for your equation in this way? Is it taken from any reference? Feb 8 '19 at 5:02
• It does seem clunky, but this is the reference: wolfram.com/language/11/partial-differential-equations/… Feb 8 '19 at 13:43
• Well, though that equation is also called wave equation, it's apparently a different one, and I won't feel surprised if the b.c. approximating infinity for that equation doesn't work on yours. In other words, there's no general way to define b.c. at infinity, AFAIK. As to the proper b.c.at infinity for your equation, I'm afraid it's off-topic here and should be asked in e.g. scicomp.stackexchange.com Feb 8 '19 at 14:49
• Somewhat related: mathematica.stackexchange.com/q/145478/1871 Feb 8 '19 at 14:50
• Thank you for the link. Your solution to mathematica.stackexchange.com/questions/128516 was also helpful and helped me see that my question here is confusing: Solutions to the above PDE are not simple waves - but the question was meant to contain the full PDE above only inside the cavity, with the simple wave equation on the outside. Perhaps I should edit the question to clarify this. From your answer in the present link, I was able to write a solution to the clarified problem which I could post here. Feb 8 '19 at 18:10

We can obtain a solution by modifying the above PDE. To do so, we introduce a smoothing function sm[x] which effectively removes the problematic functions outside of the potential. Then, thanks to @xzczd 's suggestions and the links in the comments, we implement the absorbing boundary condition (abc) for the simple wave equation.

Geometry:

{xin, xmax, tmax, pts} = {-60, 100, 300, 350};
{x0, xp, w, wp, v} = {50, 0, 2.5, 2, -1};
c = 0.4;
g[x_, t_] := Exp[-(((x - x0) - v*t)/w)^2]
gp[x_] := Exp[-((x - xp)/wp)^2]
sm[x_] := (-Tanh[2 wp x] + 1)/2

Wave equation and boundary conditions:

mol[n_Integer] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n}}
weqn = D[psi[x, t], {x, 2}] - D[psi[x, t], {t, 2}] -
2 c*I*D[psi[x, t], {t, 1}] sm[x] + (c^2 sm[x] - gp[x])*
psi[x, t] == 0;
ic = {psi[x, 0] == g[x, 0],
Derivative[0, 1][psi][x, 0] == Derivative[0, 1][g][x, 0]};
bcd = {psi[xin, t] == 0};
abcr = D[psi[x, t], x] + D[psi[x, t], t] == 0 /. x -> xmax;

Solution:

nsol = NDSolveValue[{weqn, ic, bcd, abcr},
psi, {x, xin, xmax}, {t, 0, tmax}, Method -> mol[pts]];