Boundary conditions at infinity for 1+2D wave equation in Mathematica 7

Question

To solve a waves equation, I need to define some boundary conditions. The wave is propagating on an infinite plane, and it's not a membrane fixed on some fixed support. I'm have difficulties in setting up the proper conditions with this working code:

Tmax = 20;

WaveEquation = 
  D[Phi[t, x, y], t, t] - D[Phi[t, x, y], x, x] - D[Phi[t, x, y], y, y] == 0;

PHI0[x_, y_] = Tanh[x] - 1.5Exp[-0.25(x^2 + y^2)]; (* Initial field *)

BoundaryConditions = { (* The issue is here ! *)
   Phi[t, -10, y] == PHI0[-10, y],
   Phi[t, 10, y] == PHI0[10, y],
   Phi[t, x, -10] == PHI0[x, -10],
   Phi[t, x, 10] == PHI0[x, 10]
   };

InitialConditions = {
   Phi[0, x, y] == PHI0[x, y],
   (D[Phi[t, x, y], t]/.t -> 0) == 0
   };

WaveSolution = NDSolve[
   Flatten@{WaveEquation, InitialConditions, BoundaryConditions},
    Phi,
   {t, 0, Tmax},
   {x, -10, 10},
   {y, -10, 10},
   PrecisionGoal -> 2
   ];

Manipulate[
 Plot3D[
  Evaluate[Phi[t, x, y]/.WaveSolution/.t -> time],
  {x, -10, 10},
  {y, -10, 10},
  PlotPoints -> {20, 20},
  PlotRange -> {{-10, 10}, {-10, 10}, {-3, 3}},
  Axes -> True,
  MeshFunctions -> (#3&),
  ColorFunction -> "Rainbow",
  ImageSize -> 700,
  Method -> {"RotationControl" -> "Globe"},
  SphericalRegion -> True
  ],
 {time, 0, Tmax, 0.01}
 ]

By running this code, we could clearly see the waves reflections on the bounding box. I want to get rid of that effect, since the wave is supposed to be propagating on an infinite plane. What can I do here? Take note that I'm working on Mathematica 7.0, so the option PeriodicBoundaryCondition isn't recognized by my version, and it would be inapropriate anyway since the waves would get back to the origin after a while.

So how can I tell Mathematica that there isn't any boundary to constraint the waves?

---

EDIT1: Is there a way to tell Mathematica that all waves should get truncated at the boundaries, without producing any reflections, creating something equivalent to a boundary at infinity? My problem mainly comes from the reflections at the borders {{-10, 10}, {-10, 10}} that I would like to prevent.

---

EDIT2: My question is similar to this one (but without an answer yet): Absorbing/Derivative Boundary Conditions

Answer 1

Aha, finally here comes a chance to share my implementation of 2nd order absorbing boundary condition (ABC). To understand why the b.c. is implemented in this way, one may refer to Chapter 6 of Lloyd N. Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations but to be honest I myself don't understand the underlying math very well, either.

I've used pdetoode to discretize the PDE to a system of ODE because NDSolve cannot directly handle 2nd order ABC.

The code is tested in v6.0.0 so should work in v7, too.

domain = {-10, 10}; tend = 20;

patt = x | y | t;
Table[{Numerator@
    Together@With[{theothervar = First@Complement[{x, y}, {var}]}, 
      var - PadeApproximant[#, {s, 0, {2, 0}}] & /@ 
        Simplify[var /. Solve[t^2 == x^2 + y^2, var] /. theothervar -> s t, t > 0] /. 
       s -> theothervar/t], var == Reverse@domain // Thread}\[Transpose], {var, {x, y}}]

bc = % /. {(var1 : patt) (var2 : patt) :> D[u[t, x, y], var1, var2], 
           (var : patt)^n_ :> D[u[t, x, y], {var, n}]}

With[{u = u[t, x, y]}, 
 ic = {D[u, t] == 0, u == Tanh[x] - 1.5 Exp[-0.25 (x^2 + y^2)]} /. t -> 0;
 eq = Piecewise[Flatten[bc, 1], D[u, t, t] - (D[u, x, x] + D[u, y, y])] == 0]

points = 50; difforder = 2;
grid = Rescale[Range@points, {1, points}, domain];

(*Definition of pdetoode isn't included in this post,
  please find it in the link above.*)
ptoofunc = pdetoode[u[t, x, y], t, {grid, grid}, difforder];
ode = ptoofunc@Simplify`PWToUnitStep@eq;
odeic = ptoofunc@ic;

sollst = Partition[
   NDSolve[{ode, odeic}, Flatten@Outer[u, grid, grid], {t, 0, tend}][[1, All, -1]], 
   points];

solfunc = rebuild[sollst, {grid, grid}]

Table[Plot3D[solfunc[t, x, y], {x, ##}, {y, ##}, PlotRange -> {-1.5, 1.5}] & @@ 
   domain, {t, 0, 20, 1/2}] // ListAnimate

P.S. For those in higher versions (At least $VersionNumber>=10 I guess? ), the implementation of ABC and PML in the tutorial Acoustics in the Time Domain may be a better choice.