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

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!

• In this answer, I just used a narrow Gaussian instead of a point. Maybe that's good enough...? After all, it's all discretized anyway, so points are never points. – Jens Jan 16 '15 at 23:57
• @Jens Yeah, that sounds reasonable. Maybe you can consider giving an answer? – xzczd Jan 17 '15 at 1:25

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["NDSolveFEM"];
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.

• I think you're right about the continuous limit being ill-defined, and that's why I modified the problem to put the driving in as a source term of finite size. The question still makes sense if we consider the grid as a given fixed structure, but this means NDSolve by itself doesn't seem to be the right tool here. – Jens Jan 17 '15 at 19:25
• Interesting……Well, does this mean if we stick a infinite thin rod on an unbroken membrane, the membrane won't have any displacement except on the stuck point? – xzczd Jan 19 '15 at 3:17
• @xzczd: Well, the Laplace equation approximates the displacement of an elastic membrane only when the gradient of the displacement is small. What you could argue is that in the limit as the displacement of the rod goes to zero, so does the relative size of the region it affects. The Laplace equation is a better model for heat distribution: if you stick an infinitely thin pin held at a constant warm temperature on a colder membrane, it won't heat up the membrane at all. – user484 Jan 19 '15 at 4:25
• Concerning your question if a single point can be used as a boundary condition: Yes, if the condition is a Dirichlet condition and not if the condition is a Neumann condition. Dirichlet conditions act on points (in the mesh) and Neumann conditions act on edges over the region. Good answer. If you have suggestions for the FEM documentation let me know. – user21 Jan 19 '15 at 10:13
• @Rahul, even though "IncludePoints" is undocumented (you can hang/crash your system with it) I fixed this in the development version. The include point is now present in the PoineElement list. – user21 Jan 28 '15 at 16:09

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.