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.

  • 2
    $\begingroup$ An eigendecomposition of the finite element discretization of the Laplacian should give you what you're looking for $\endgroup$ Commented Jul 7, 2022 at 19:09
  • $\begingroup$ Do you know of any references which would discuss this? $\endgroup$
    – vibe
    Commented Jul 7, 2022 at 20:05
  • $\begingroup$ If we want to "estimate" those eigenfunctions rather than solve for them, I would consider a formulation of the Laplace equation as a minimization problem, and try to find an approximate solution there it in the style of the variational approach of estimating the ground state energy in quantum mechanics. $\endgroup$ Commented Jul 8, 2022 at 5:36
  • $\begingroup$ @Maxim: I don't necessarily want to solve the Laplace equation. I want a method which will produce all of the eigenfunctions which can then be used to construct a full solution later. $\endgroup$
    – vibe
    Commented Jul 8, 2022 at 16:58
  • $\begingroup$ I was meaning solving for the eigenfunctions (although with the knowledge of eigenfunctions you effectively have solved the Laplacian equation as well). If the question is about numerical solution for the eigenfunctions, what is needed to do there is discretize the Laplace operator with your favorite method (finite element or finite difference, does not matter), and then the operator becomes a matrix. Then you put this matrix in a standard package for finding eigenvalues and vectors for matrices, that's it. $\endgroup$ Commented Jul 9, 2022 at 4:21

2 Answers 2


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.

  • $\begingroup$ Thank you for your comment. Do you know of any references, example code, etc which show how to actually perform the steps in Theorem 3? $\endgroup$
    – vibe
    Commented Jul 8, 2022 at 19:01
  • $\begingroup$ @vibe, there are several answers within this site for that. See for example this one $\endgroup$
    – nicoguaro
    Commented Jul 8, 2022 at 22:41
  • $\begingroup$ That's right -- the eigenfunctions are simply the functions associated with the eigenvectors of the discretized operator. $\endgroup$ Commented Jul 9, 2022 at 1:52
  • $\begingroup$ In order for the $A_i$ to be non-trivial, it seems required that $f \neq 0$. So does this mean that this method will not work for the Laplace equation, where $f = 0$? $\endgroup$
    – vibe
    Commented Jul 9, 2022 at 16:29
  • $\begingroup$ No, there are also terms that relate to the boundary values $g$. I updated the answer accordingly. $\endgroup$ Commented Jul 10, 2022 at 23:16

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.