Least effort to handle a point source inside the domain of PDE(s)

Question

By point source I mean a constrained condition at one point inside the domain of PDE(s). For example:

$$\frac{\partial ^2u(t,x,y)}{\partial t^2}=\frac{\partial ^2u(t,x,y)}{\partial x^2}+\frac{\partial ^2u(t,x,y)}{\partial y^2}$$ $$u(t,0,0)=\sin (10 t)$$ $$u(0,x,y)=0,u^{(1,0,0)}(0,x,y)=0$$ $$u(t,-1,y)=0,u(t,1,y)=0$$ $$u(t,x,-1)=0,u(t,x,1)=0$$ $$0\leq t\leq 3,-1\leq x\leq 1,-1\leq y\leq 1$$

This can model… Er… a square edge-fixed membrane with one tip of an ultra thin moving rod stuck on the center. The condition $u(t,0,0)=\sin (10 t)$ is exactly a point source. (Notice it's not completely compatible with the initial conditions but it's not a big deal. )

NDSolve can't solve this problem directly (at least now):

NDSolve[{
  D[u[t, x, y], t, t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y],
  u[t, 0, 0] == Sin[10 t],
  u[0, x, y] == 0,
  Derivative[1, 0, 0][u][0, x, y] == 0,
  u[t, -1, y] == 0,
  u[t, 1, y] == 0,
  u[t, x, -1] == 0,
  u[t, x, 1] == 0},
 u, {t, 0, 3}, {x, -1, 1}, {y, -1, 1}]

NDSolve::bcedge: Boundary condition u[t,0,0]==Sin[t] is not specified on a single edge of the boundary of the computational domain. >>

Of course FDM can handle the point source naturally:

ans = (Reap@
      With[{n = 50, c = 1}, 
        With[{dx = (1 - (-1))/(n - 1), Courant = Sqrt[2]/2}, 
         With[{dt = (Courant dx)/c}, 
          Compile[{}, 
           Module[{z1, z2}, z1 = z2 = Table[0., {n}, {n}]; 
            Do[{z1, z2} = {z2, z1}; 
             z1[[Ceiling[n/2], Ceiling[n/2]]] = Sin[10 t]; 
             Do[z2[[i, j]] = 
               z1[[i, j]] + z1[[i, j]] - z2[[i, j]] + 
                Courant^2 (z1[[i - 1, j]] + z1[[i + 1, j]] + 
                   z1[[i, j - 1]] + z1[[i, j + 1]] - 
                   4 z1[[i, j]]), {i, 2, n - 1}, {j, 2, n - 1}]; 
             Sow[z1], {t, 0, 3, dt}]]]]]][])[[-1, 1]];
ListPlot3D[#, Mesh -> False, PlotRange -> {-1, 1}] & /@ ans;
SystemOpen@Export["a.gif", %];

But can we more or less benefit from NDSolve or other existed tools in Mathematica, instead of doing something from scratch? Is NDSolve completely useless in this situation?

A general solution is the best, but opportunistic ones, I mean, solutions that are only suited for the specific example above are also welcomed!

Answer 1

The good news is that yes, there is an easy way to put your problem into NDSolve by using the new finite element functionality in v10. The bad news is that it seems the specific problem you're trying to solve is ill-posed.

---

NDSolve can now handle internal boundaries; see e.g. the first figure under "Details" for DirichletCondition. Generating a mesh with such internal boundaries is described in the "Element Mesh Generation" tutorial. I don't know if a single constrained point technically counts as a "boundary", but it seems to work.

Create a spatial mesh with a node at the point source:

Needs["NDSolve`FEM`"];
bmesh = ToBoundaryMesh[
   "Coordinates" -> {{-1, -1}, {-1, 1}, {1, 1}, {1, -1}, {0, 0}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}]}];
mesh = ToElementMesh[bmesh];
Show[mesh["Wireframe"], Graphics[{Red, PointSize[Large], Point[{0, 0}]}]]

(One could also use the not-really-documented "IncludePoints" option, as in this other answer.)

Direcly specifying u[t, 0, 0] doesn't work, as you already know, but DirichletCondition does:

sol = NDSolve[{
   D[u[t, x, y], t, t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y], 
   DirichletCondition[u[t, x, y] == Sin[10 t], x == 0 && y == 0], 
   u[0, x, y] == 0,
   Derivative[1, 0, 0][u][0, x, y] == 0, 
   u[t, -1, y] == 0,
   u[t, 1, y] == 0,
   u[t, x, -1] == 0, 
   u[t, x, 1] == 0},
  u, {t, 0, 3}, {x, y} ∈ mesh];

It complains that "NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified", which is to be expected. But it gives a solution anyway.

frames = Table[
   Plot3D[u[t, x, y] /. sol, {x, -1, 1}, {y, -1, 1}, 
    PlotRange -> {-1, 1}, Mesh -> None, PlotStyle -> White], {t, 0, 3, 0.05}];
Export["a.gif", frames];

---

We start to see a problem if we change the resolution of the mesh. To avoid the massive memory requirements of a uniformly refined mesh, it's better to refine only where the solution changes rapidly, i.e. in the neighbourhood of the point source. One can use

mesh = ToElementMesh[bmesh, 
   "MeshRefinementFunction" -> Function[{vertices, area}, 
     area > Max[a, 1*^-2 Min[Norm /@ vertices]]]];

which smoothly refines the mesh to have elements of area $a$ near the point source at the origin. Here are some meshes with $a=10^{-2}$, $10^{-4}$, and $10^{-6}$, followed by zooms to $[-0.1,0.1]\times[-0.1,0.1]$:

And here are the corresponding solutions at $t=1$:

The solutions seem to be getting weaker the finer we make the mesh. What's going on? I don't know for absolutely certain, but I'm guessing that the problem is ill-posed and the solution we've computed is essentially an artifact of the numerical discretization. As an analogy, consider the Laplace problem on a punctured domain with Dirichlet boundary conditions: $$\begin{align} \nabla^2f(x,y)&=0&\text{for }&x\in\Omega\setminus\{(0,0)\},\\ f(x,y)&=0&\text{for }&x\in\partial\Omega,\\ f(0,0)&=1. \end{align}$$ You can solve this numerically and obtain a reasonable-looking numerical solution, but it is an illusion because a one-point set has zero capacity for the Laplacian, and if you refine the mesh the solution goes to zero. I believe the same thing is happening here. Numerically, the energy that the source imparts to the system is mesh-dependent, being related to the area of its neighbouring elements. Theoretically, I guess there is no solution.

---

So yeah. Can you use NDSolve for this problem? You can. But... maybe you shouldn't.

Disclaimer: I am not a functional analyst and this is not mathematical advice. Consult your friendly neighbourhood applied mathematician.

Answer 2

The periodic driving at one point doesn't seem to be compatible with the boundary conditions expected by NDSolve, so I modified the problem in two ways: first, broaden the point source into a Gaussian, and then incorporate this driving as a source term in the actual differential equation.

So we're actually solving the inhomogeneous wave equation here. For the initial condition, I chose something that looks like the driving term spatially. The result is still not very satisfactory because NDSolve is much slower than your hand-coded finite-difference scheme.

σ = .1;

sol = u /. First@NDSolve[{
     D[u[t, x, y], t, t] == 
      D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + 
       Sin[10 t] Exp[-(x^2 + y^2)/(2 σ^2)],
     Derivative[1, 0, 0][u][0, x, y] == 0,
     u[0, x, y] == .004 Exp[-(x^2 + y^2)/(2 σ^2)],
     u[t, -1, y] == 0,
     u[t, 1, y] == 0,
     u[t, x, -1] == 0,
     u[t, x, 1] == 0},
    u,
    {t, 0, 3},
    {x, -1, 1},
    {y, -1, 1}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 100}}];

frames = 
  Table[Plot3D[sol[t, x, y], {x, -1, 1}, {y, -1, 1}, 
    PlotRange -> .02 {-1, 1}, PlotPoints -> 20, Mesh -> False], {t, 0,
     3, .05}];

ListAnimate[frames]

The parameter $\sigma$ dictates the spatial width of the Gaussian in the driving term, and I added the "Method" options explicitly so that we can control the value of "MaxPoints" if you decide to make $\sigma$ smaller than what I chose. The smaller $\sigma$, the larger "MaxPoints" needs to be. This makes NDSolve painfully slow if you really want to approach a well-localized point source. Of course, the plotting of the resulting InterpolatingFunction eats up some time, too. So in terms of speed, your discretized calculation is faster both in the actual differential-equation solving and in the plotting steps.

Nevertheless, the above at least shows how we can make NDSolve produce the desired result in principle.