Numerical estimation of eigenfunctions of Laplacian

Question

Consider the Laplace equation, $$ \nabla^2 f(r,\theta,\phi) = 0 $$ in spherical coordinates. We know that the solution to this equation can be derived analytically, and is given by, $$ f(r,\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^l \left( A_l^m \frac{1}{r^{l+1}} + B_l^m r^l \right) Y_l^m(\theta,\phi) $$ where the coefficients $A_l^m,B_l^m$ would be determined by boundary conditions on the solution $f$. This solution is built out of the eigenfunctions of the Laplacian $\frac{1}{r^{l+1}} Y_l^m$ and $r^l Y_l^m$.

My question is, do numerical methods exist which can generate these basic building blocks $\frac{1}{r^{l+1}} Y_l^m$ and $r^l Y_l^m$ if we didn't already know them from analytic results? Ultimately, I want to build these types of functions for a coordinate system where the Laplacian does not separate and for which there are no known analytic solutions.

Answer 1

Fundamentally, here are the building blocks of what you are asking for:

• Consider solving the problem $$ -\Delta u = f $$ in a domain $\Omega$ with boundary values $u=g$ on $\partial\Omega$.

• Theorem 1: For any $\Omega,f,g$, the solution can be written as $$ u(\mathbf x) = \sum_{i=1}^\infty A_i \varphi_i(\mathbf x) $$ where the functions $\varphi_i$ are the eigenfunctions of the Laplace operator and satisfy $$ -\Delta \varphi_i = \lambda_i \varphi_i, $$ and the $A_i$ are expansion coefficients that can be computed via $A_i = \frac{1}{\lambda_i}\int_\Omega \varphi_i(\mathbf x) f(\mathbf x)\, dx$ plus some terms that come from the boundary values $g$. (The theorem is true because the eigenfunctions of the Laplace operator form a complete basis of $L_2(\Omega)$ and consequently of the solution space in which $u$ lies.)

• Theorem 2: For a given $\Omega,f,g$, the coefficients $A_i$ must decay, that is: $A_i\rightarrow 0$. As a consequence, you can approximate the solution by truncating the sum as $$ u(\mathbf x) \approx \sum_{i=1}^N A_i \varphi_i(\mathbf x). $$ (This theorem is true because the norm of $u$ must be finite, and that can only be the case if the sum $\sum_{i=1}^\infty A_i$ converges -- which it only does if the coefficients decay sufficiently fast.)

• Theorem 3: The eigenfunctions $\varphi_i$ can be approximated numerically, for example using the finite element method, by functions $\varphi_i^h$. Then the solution of the original problem can be approximated as $$ u(\mathbf x) \approx u^h(\mathbf x) = \sum_{i=1}^N A_i \varphi_i^h(\mathbf x). $$

I will note, however, that using this expansion into eigenfunctions is not an efficient way to find an approximate solution of the Laplace equation.

Answer 2

The solution of $$ \nabla^2 f(r,\theta,\phi) = 0 $$ can be written in terms of eigenfunctions of the Laplacian. As you said, the expansion coefficients then should be chosen as to satisfy the differential equation as well as the boundary conditions (this works as the eigenfunctions form a complete set of the Hilbert space).

In spherical coordinates, the expansion into eigenfunctions is given by $$ f(r,\theta,\phi) = \sum_{klm} R_{lm}(r) Y_l^m(\theta,\phi) $$ The functions $R_{lm}$ satisfy the radial Laplace equation, $$ \frac{1}{r^2}\frac{d}{dr} \left(r^2\frac{dR_{lm}(r)}{dr}\right) - \frac{\ell(\ell+1)}{r^2} R_{lm}(r)= E_{lm} R_{lm}(r) $$ The solutions you mentioned are eigenfunctions of this equation. For the Laplace equation, the solution actually does not depend on $m$ due to rotational symmetry around the $z$-axis (so you could also use an expansio into Legendre polynomials than into sperical harmonics, but the latter is more general and also applicable to non-symmetric problems).

As mentioned by @WolfgangBangerth, one is looking for the regular component of the solution to the above equation (which is the one with the $A$-coefficients). For this, you can make another ansatz by using $u_{lm}(r) =r R_{lm}(r)$, which leads to the equation

$$ -\frac{d^2u_{lm}(r)}{dr^2} + \frac{\ell(\ell+1)}{r^2} u_{lm}(r)= E_{lm} u_{lm}(r) $$

It is this one-dimensional equation which you should solve numerically in order to arrive at the spherical eigenfunctions of the Laplacian. You can do that by the usual methods, e.g. finite differences, -elements, or spectral methods.

Further, by looking at the form of the equation, it directly becomes clear that the (inverse) polynomial functions in $r$ you mentioned are in fact eigenfunctions.