You can't solve Laplace's equation with boundary conditions on isolated points. But why?

Question

Consider a bounded region $\Omega\subset\mathbb R^n$ with a finite number of "holes" $X=\{x_1,\ldots,x_k\}$ that are isolated points in its interior. I'm pretty sure that in more than one dimension, it doesn't make sense to solve Laplace's equation with Dirichlet boundary conditions on $X$. That is, there does not exist any "valid"* $f$ such that $$\begin{align} \nabla^2f(x) &= 0 & \text{for } & x\in\operatorname{int}\Omega\setminus X, \\ f(x) &= 0 & \text{for } & x \in \partial\Omega, \\ f(x_i) &= y_i \ne 0 & \text{for } & i=1,\dots,k. \end{align}$$ I feel this is related to the fact that the Green's function for the Laplacian has a singularity at the origin, but I don't know how to get from here to there. Is there a straightforward proof, or a well-known theorem that the desired result follows from directly?†

Also, it would be useful to know if there are generalizations to higher-order differential equations like the biharmonic equation, $\nabla^4f=0$.

Update: As the comments on a now-deleted answer imply, this boils down to the fact that any not-too-singular point discontinuity in a harmonic function is removable. Wikipedia states this fact without proof or references, but mentions that Riemann's theorem proves the analogous theorem for holomorphic functions on the complex plane. What's the proof for harmonic functions on $\mathbb R^{n\ge2}$?

---

*I say "valid" because you could have a function that's identically zero everywhere except on $X$, but that "shouldn't count" because you could satisfy absolutely any Dirichlet boundary condition that way. Please let me know the correct PDE-theoretic term for this.

†I've had this question come up multiple times, whenever someone thinks of doing scattered data interpolation by just plugging the data points into Laplace's equation. Numerically, they end up with a spiky-looking function; mathematically, the solution doesn't even exist. In the future, I would like to be armed with a theorem when I have to explain to them why this cannot possibly work.

Answer 1

Reference

I recommend Classical Potential Theory by Armitage and Gardiner, section 5.2. Corollary 5.2.3 implies that if a bounded function is harmonic in $\Omega\setminus F$ where $F$ is a finite set, then it extends as a harmonic function in $\Omega$. In your setup, the extended function would contradict the maximum principle.

Proof

I will use Green's function, as you expected. If $n>2$, form the sum $$ f(x) - \epsilon\sum_{i=1}^k \|x-x_k\|^{2-n} $$ Note that $f$ tends to $-\infty$ at the points $x_i$. Letting $f(x_i)=-\infty$, we obtain a subharmonic function in all of $\Omega$. Indeed, one characterization of subharmonic functions is the local sub-mean value property: for every point $a$ there is $r_a>0$ such that the average of $f$ on every sphere $S(a,r)$ with $r<r_a$ is at least $f(a)$. This property holds at $x_i$ trivially, and at other points by the harmonicity of $f$.

Subharmonic functions satisfy the maximum principle. For a fixed $p\in \Omega\setminus \{x_i\}$ it says
$$ f(p) \le O(\epsilon) $$ and since $\epsilon$ was arbitrarily small, $f\le 0$ in $\Omega$. The same consideration applies to $-f$; thus, $f\equiv 0$.

For $n=2$ the proof is the same, but has the logarithm in it: $$ f(x) + \epsilon\sum_{i=1}^k \log\|x-x_k\| $$

Other PDE

The keyword is capacity: does a one-point set have zero capacity for the given differential operator? The Laplace equation is the Euler-Lagrange equation for the functional $\int|\nabla f|^2$, so the relevant capacity of point $p$ is $$ \inf \left\{\int|\nabla f|^2 : f=0\text{ on boundary, } f(p)=1 \right\} $$ The infimum can be taken over all smooth functions, or over all Lipschitz functions: does not matter. The capacity of a point is zero, as you can see by letting $f$ be a scaled and truncated version of Green's function. Since the solution of boundary value problem is meant to minimize the functional for the given boundary data, the fact that the infimum is zero says there isn't a solution.

For the biharmonic equation, the relevant functional is $\int|\Delta f|^2$. Is it true that $$ \inf \left\{\int|\Delta f|^2 : f=0\text{ on boundary, } f(p)=1 \right\} =0 \quad ? $$ Yes in dimensions $n\ge 4$, because biharmonic Green's function is unbounded there (source). No in dimensions $n\le 3$. The reason is that the $L^2$ norm of the Laplacian controls $W^{2,2}$ norm (roughly speaking), and in dimensions $n\le 3$ the Sobolev embedding implies that $W^{2,2}$ norm controls pointwise values of a function.

So, the biharmonic equation can be used for interpolation in low dimensions. In two dimensions, this is known as a thin plate spline. More generally, there are polyharmonic splines: the higher the dimension, the more one iterates the Laplacian to keep Green's function bounded.