# Preface

This section presents the basic properties of solutions to the Laplace and Helmholtz equations.

Introduction to Linear Algebra with Mathematica

# Laguerre Expansion

The Laguerre polynomials, named after a French mathematician Edmond Nicolas Laguerre (1834--1886), are solutions of Laguerre's equation

$x\,y'' + \left( 1-x \right) y' +n\,y =0, \quad\mbox{or} \quad \frac{{\text d}}{{\text d} x} \left( x\,e^{-x} \,y' \right) + n\,e^{-x}\,y =0 , \qquad x\in (0,\infty ).$
This second order linear differential equation and its solutions were first analyzed by Pafnuty Chebyshev (1821--1894) in 1859 almost 20 years before Laguerre studied them. This equation has a polynomial solution, which may be defined by the Rodrigues formula (the French banker Olinde Rodrigues (1795--1851) defined it in 1816)
$L_n (x) = \frac{e^x}{n!} \, \frac{{\text d}}{{\text d} x} \left( x^n e^{-x} \right) , \qquad n=0,1,2,\ldots .$
These polynomials are called Chebyshev--Laguerre or simply Laguerre polynomials and they can be written in closed form:
$L_n (x) = \sum_{k=0}^n \binom{n}{k} \, \frac{(-1)^k}{k!} \, x^k .$
Generating function:
$\frac{1}{1-t} \, e^{-xt/(1-t)} = \sum_{n\ge 0}^n L_n (x)\, t^n .$
One can also define the Laguerre polynomials recursively,
$L_{n+1} (x) = \frac{1}{n+1} \left[ \left( 2n+1-x \right) L_n (x) -n\,L_{n-1} \right] , \quad L_0 (x) =1, \quad L_1 (x) = 1-x, \quad n=2,3,\ldots .$
Arbitrary function for which the integral
$\| f \|^2 = \int_0^{\infty} e^{-x} |f(x)|^2 {\text d} x < \infty$
is finite, can be extended into Laguerre series:
$f(x) = \sum_{i\ge 0} f_i L_i (x) , \qquad f_i = \int_0^{\infty} L_i (x) \,e^{-x} f(x) \, {\text d} x ,$
which converges in $$L^p (-\infty , \infty )$$ for $$p \in \left( \frac{4}{3} , 4 \right) .$$ Such expansion is based on the orthogonal property of Laguerre polynomials:
$\int_0^{\infty} L_n (x) \,L_m (x) \, e^{-x} \, {\text d} x = \delta_{n,m} \equiv \begin{cases} 1 , & \quad \mbox{if n=m}, \\ 0 , & \quad \mbox{if n\ne m}. \end{cases}$

More generally, the associated Laguerre equation

$x\,y'' + \left( \alpha + 1-x \right) y' +n\,y =0, \quad\mbox{or} \quad \frac{{\text d}}{{\text d} x} \left( x^{1+\alpha}\,e^{-x} \,y' \right) + n\,x^{\alpha} \,e^{-x}\,y =0 , \qquad x\in (0,\infty ).$
has polynomial solutions, denoted by $$L_n^{(\alpha )} (x)$$ and called the generalized/associated Laguerre polynomials, or Sonine polynomials, after their inventor, a Russian mathematician Nikolay Yakovlevich Sonin (1849--1915). They can be defined either by the Rodrigues formula
$L_n^{(\alpha )} (x) = x^{-\alpha} \,\frac{e^x}{n!} \, \frac{{\text d}}{{\text d} x} \left( x^{n+\alpha} e^{-x} \right) , \qquad n=0,1,2,\ldots .$
or through its ordinary generating function
$\frac{1}{(1-t)^{\alpha +1}} \, e^{-xt/(1-t)} = \sum_{n\ge 0}^n L_n^{(\alpha )} (x)\, t^n .$

Example: Consider a power function $$f(x) = 4\,x^3 -1 .$$ This function has a finite square norm with weight $$e^{-x} :$$

$\| 4\,x^3 -1 \|_2^2 = \int_0^{\infty} \left( 4\,x^3 -1 \right)^2 \, e^{-x} \, {\text d} x = 11473.$
Therefore, this function can be expanded into convergent Laguerre series (which is actually a finite sum):
$4\,x^3 -1 = \sum_{k\ge 0} c_k L_k (x) ,$
where
$c_0 = \int_0^{\infty} \left( 4\,x^3 -1 \right) \,e^{-x} \, {\text d} x = 23 , \quad c_1 = -72, \quad c_2 =72, \quad c_3 = -24 .$
All other coefficients are zeroes, and we get the identity:
$4\,x^3 -1 = 23\,L_0 (x) -72\,L_1 (x) + 72\,L_2 (x) -24\,L_3 (x) ,$

Example: Consider a power function $$\displaystyle f(x) = \frac{4\,x^3 -1}{x^2 +1} .$$ This function has a finite square norm with weight $$e^{-x} :$$

$\left\| \frac{4\,x^3 -1}{x^2 +1} \right\|_2^2 = \int_0^{\infty} \left( \frac{4\,x^3 -1}{x^2 +1} \right)^2 \, e^{-x} \, {\text d} x \approx 21.3606.$
Therefore, this function can be expanded into convergent Laguerre series
$\frac{4\,x^3 -1}{x^2 +1} = \sum_{k\ge 0} c_k L_k (x) ,$
where
\begin{align*} c_0 &= \int_0^{\infty} \left( \frac{4\,x^3 -1}{x^2 +1} \right) \, e^{-x} \, {\text d} x \approx 2.00504, \\ c_1 &= \int_0^{\infty} \left( \frac{4\,x^3 -1}{x^2 +1} \right) \, e^{-x} \,L_1 (x) \, {\text d} x \approx -4.13738, \\ c_2 &= \int_0^{\infty} \left( \frac{4\,x^3 -1}{x^2 +1} \right) \, e^{-x} \,L_2 (x) \, {\text d} x \approx 0.217678, \\ c_3 &= \int_0^{\infty} \left( \frac{4\,x^3 -1}{x^2 +1} \right) \, e^{-x} \,L_3 (x) \, {\text d} x \approx 0.260623, \end{align*}
and so on, getting $$c_4 \approx 0.223731, \ c_5 \approx 0.17164.$$ Now we build Laguerre approximation with six terms:
c3 = NIntegrate[ LaguerreL[3, x]*(4*x^3 - 1)*Exp[-x]/(x*x + 1), {x, 0, Infinity}]
laguerre = c0 + c1*LaguerreL[1, x] + c2*LaguerreL[2, x] + c3*LaguerreL[3, x] + c4*LaguerreL[4, x] + c5*LaguerreL[5, x]
Plot[{(4*x^3 - 1)/(x*x + 1), laguerre}, {x, 0, 5}, PlotStyle -> {Blue, Orange}]

Example: Find the Fourier-Laguerre series expansion of the power function $$f(x) = x^p , \ p>1 .$$

To answer this question, we need to find coefficients ck in the Laguerre expansion:

$c_k = \int_0^{\infty} x^p e^{-x} L_k (x)\,{\text d} x = \frac{1}{k!} \, \int_0^{\infty} x^p \,\frac{{\text d}^k}{{\text d} x^k} \left( x^k e^{-x} \right) , \quad k=0,1,2,\ldots .$
Starting with k = 0, we have
$c_0 = \int_0^{\infty} x^p e^{-x} \,{\text d} x = \Gamma (p+1) ,$
where Γ is the gamma function of Euler. For k > 0, we integrate by parts in the integral
\begin{align*} c_k &= \frac{1}{k!} \, \int_0^{\infty} x^p \,\frac{{\text d}^k}{{\text d} x^k} \left( x^k e^{-x} \right) \,{\text d} x = \left. \frac{1}{k!} \, x^p \,\frac{{\text d}^{k-1}}{{\text d} x^{k-1}} \left( x^k e^{-x} \right) \right\vert_{x=0}^{\infty} - \frac{p}{k!} \, \int_0^{\infty} x^{p-1} \,\frac{{\text d}^{k-1}}{{\text d} x^{k-1}} \left( x^k e^{-x} \right) \,{\text d} x \\ &= -\left. \frac{p}{k!} \, x^{p-1} \,\frac{{\text d}^{k-2}}{{\text d} x^{k-2}} \left( x^k e^{-x} \right) \right\vert_{x=0}^{\infty} + \frac{p(p-1)}{k!} \, \int_0^{\infty} x^{p-2} \,\frac{{\text d}^{k-2}}{{\text d} x^{k-2}} \left( x^k e^{-x} \right) \,{\text d} x \\ &= (-1)^k \,\frac{1}{k!} \, p^{\underline{k}} \, \Gamma (p+1) , \end{align*}
where $$p^{\underline{k}} = p(p-1) \cdots (p-k+1)$$ is kth falling factorial. Hence, the Fourier-Laguerre series expansion of the power function is given by
$x^p = \Gamma (p+1) + \Gamma (p+1) \,\sum_{k\ge 1} \frac{(-1)^k}{k!} \, p^{\underline{k}} \, L_k (x) .$
If p is a positive integer, the above series becomes a polynomial of degree p because falling factorial $$p^{\underline{k}} =0$$ for k > p. Also $$\Gamma (p+1) = p!$$ for positive integer p.

Example: Dirac delta function

$\delta (x-a) = e^{-(x+a)/2} \,\sum_{k\ge 0} \, L_k (x) \, L_k (a) .$

Example: The most important application of the Laguerre polynomials is in the solution of the Schrödinger equation for the hydrogen atom. This equation is

$-\frac{\hbar^2}{2m}\, \nabla^2 \psi - \frac{Ze^2}{4\pi \epsilon_0 r} \, \psi = E\,\psi ,$
in which Z = 1 for hydrogen, 2 for ionized helium, and so on. Separating variables, we find that the angular dependence of ψ is the spherical harmonic $$Y_L^M \left( \theta , \varphi \right) .$$ The radial part, R(r), satisfies the equation
$-\frac{\hbar^2}{2m}\, \frac{1}{r^2} \, \frac{\text d}{{\text d}r} \left( r^2 \frac{{\text d}R}{{\text d}r} \right) - \frac{Ze^2}{4\pi \epsilon_0 r} \, R + \frac{\hbar^2}{2m}\, \frac{L(L+1)}{r^2} \,R = E\,R .$
For bound states, $$R \mapsto 0 \ \mbox{as } \ r \mapsto \infty ,$$ and R is finite at the origin r = 0. To solve this Sturm--Liouville problem, we use the abbreviations
$\rho = \alpha r, \qquad \alpha^2 = - \frac{8mE}{\hbar^2}, \qquad \lambda = \frac{mZe^2}{2\pi \epsilon_0 \alpha \hbar^2}, \quad E<0.$
Then for $$\chi (\rho ) = R(\rho /\alpha ) ,$$ we get
$\frac{1}{\rho^2} \, \frac{\text d}{{\text d}\rho} \left( \rho^2 \frac{{\text d}\chi (\rho)}{{\text d}\rho} \right) + \left( \frac{\lambda}{\rho} - \frac{1}{4} - \frac{L(L+1)}{\rho^2} \right) \chi (\rho ) =0.$
Its solution is
$\rho \,\chi (\rho ) = e^{-\rho^2 /2}\, \rho^{L+1}\, L_{\lambda -L-1}^{(2L+1)} (\rho ) .$
We must restrict the parameter λ by requiring it to be an integer $$n, \ n=1,2,\ldots .$$ This is necessary because the Laguerre function of nonintegral n would diverge, which is unacceptable for our physical problem, in which $$\lim_{r\to\infty} R(r) =0.$$ This restriction on λ, imposed by our boundary condition, has the effect of quantifying the energy
$E_n = -\frac{Z^2 m}{2n^2 \hbar^2} \left( \frac{e^2}{4\pi \epsilon_0} \right)^2 .$
The negative sign reflects the fact that we are dealing here with bound states, corresponding to an electron that is unable to escape to infinity, where the Coulomb potential goes to zero. Using this result for En, we have
$\alpha = \frac{m e^2}{2\pi \epsilon_0 \hbar^2} \cdot \frac{Z}{n} = \frac{2Z}{n\,a_0} , \qquad \rho = \frac{2Z}{n\,a_0} \,r ,$
with
$a_0 = \frac{4\pi \epsilon_0 \hbar^2}{me^2} , \qquad\mbox{the Bohr radius} .$

Example: Consider piecewise step function

$f(x) = \begin{cases} 1 , & \ 0 < t < 1 , \\ -1 , & \ 1 < t . \end{cases}$
Module[{a}, coef = {};
Do[coef = Append[coef, ((0.73575888234) ((LaguerreL[x - 1, 1] -
LaguerreL[x, 1]) - ((0.13533528323) (LaguerreL[x - 1, 2] - LaguerreL[x, 2]))))], {x, 1, 200}]]
coef[[2]] lagsum[m_, x_] :=
Module[{a}, 0.36787944117 + Sum[N[coef[[a]]]*LaguerreL[a, x], {a, 1, m}]]
lg[x_] = Piecewise[{{1, 0 <= x < 1}, {-1, 1 < x <= 2}, {0, x > 2}}];
lagraph[m_] := Plot[{lg[x], lagsum[m, x]}, {x, 0, 2}, PlotRange -> {-1.3, 1.3},
PlotStyle -> {{RGBColor[0, 0, 1], Thickness[0.005]}, {RGBColor[1, 0, 0], Thickness[0.005]}}]
lagraph[50]
lagraph[200]
Module[{a}, coef = {};
Do[coef =
Append[coef, ((0.73575888234) ((LaguerreL[x - 1, 1] - LaguerreL[x, 1]) - ((0.13533528323) (LaguerreL[x - 1, 2] - LaguerreL[x, 2]))))], {x, 1, 200}]]
coef[[2]] lagsum[m_, x_] :=
Module[{a},
0.36787944117 + Sum[N[coef[[a]]]*LaguerreL[a, x], {a, 1, m}]] lg[x_] = Piecewise[{{1, 0 <= x < 1}, {-1, 1 < x <= 2}, {0, x > 2}}];
lagraph[m_] := Plot[{lg[x], lagsum[m, x]}, {x, 0, 2}, PlotRange -> {-1.3, 1.3}, PlotStyle -> {{RGBColor[0, 0, 1], Thickness[0.005]}, {RGBColor[1, 0, 0], Thickness[0.005]}}]
"Plot of piecewise and laguerre with n=50 terms"
lagraph[50]
"Plot of piecewise and laguerre with n=200 terms"
lagraph[200]
"Plot of piecewise-laguerre with n=50 terms"
f[x_] = lg[x] - lagsum[50, x];
Plot[f[x], {x, 0, 2}, PlotRange -> {-1.3, 1.3}, PlotStyle -> {Thick, Red}]
"Plot of piecewise-laguerre with n=200 terms"
f2[x] = lg[x] - lagsum[200, x];
Plot[f2[x], {x, 0, 2}, PlotRange -> {-1.3, 1.3}, PlotStyle -> {Thick, Red}]
"Max of difference with n=50 terms"
Maximize[{f[x], 0 <= x <= 2}, x] "Min of difference with n=50 terms" Minimize[{f[x], 0 <= x <= 2}, x]
"Max of difference with n=200 terms"
Maximize[{f2[x], 0 <= x <= 2}, x]
"Min of difference with n=200 terms"
Minimize[{f2[x], 0 <= x <= 2}, x]
"Integral of (piecewise-laguerre)^2 over interval [0,2] for n=50 terms"
Integrate[(f[x]^2), {x, 0, 2}]
"Integral of (piecewise-laguerre)^2 over interval [0,2] for n=200 terms"
Integrate[(f2[x]^2), {x, 0, 2}]
■