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
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.
Comments
Post a Comment