17.4 A Worked Example: Spheroidal Harmonics 773
Here m is an integer, c is the “oblateness parameter,” and λ is the eigenvalue. Despite
the notation, c2 can be positive or negative. For c2 > 0 the functions are called
“prolate,” while if c2 < 0 they are called “oblate.” The equation has singular points
at x = ±1 and is to be solved subject to the boundary conditions that the solution be
regular at x = ±1. Only for certain values of λ, the eigenvalues, will this be possible.
If we consider ﬁrst the spherical case, where c = 0, we recognize the differential
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
m
equation for Legendre functions Pn (x). In this case the eigenvalues are λmn =
n(n + 1), n = m, m + 1, . . . . The integer n labels successive eigenvalues for
ﬁxed m: When n = m we have the lowest eigenvalue, and the corresponding
eigenfunction has no nodes in the interval −1 < x < 1; when n = m + 1 we have
the next eigenvalue, and the eigenfunction has one node inside (−1, 1); and so on.
A similar situation holds for the general case c2 = 0. We write the eigenvalues
of (17.4.1) as λmn (c) and the eigenfunctions as Smn (x; c). For ﬁxed m, n =
m, m + 1, . . . labels the successive eigenvalues.
The computation of λmn (c) and Smn (x; c) traditionally has been quite difﬁcult.
Complicated recurrence relations, power series expansions, etc., can be found
in references [1-3]. Cheap computing makes evaluation by direct solution of the
differential equation quite feasible.
The ﬁrst step is to investigate the behavior of the solution near the singular
points x = ±1. Substituting a power series expansion of the form
∞
S = (1 ± x)α ak (1 ± x)k (17.4.2)
k=0
in equation (17.4.1), we ﬁnd that the regular solution has α = m/2. (Without loss
of generality we can take m ≥ 0 since m → −m is a symmetry of the equation.)
We get an equation that is numerically more tractable if we factor out this behavior.
Accordingly we set
S = (1 − x2 )m/2 y (17.4.3)
We then ﬁnd from (17.4.1) that y satisﬁes the equation
d2 y dy
(1 − x2 ) 2
− 2(m + 1)x + (µ − c2 x2 )y = 0 (17.4.4)
dx dx
where
µ ≡ λ − m(m + 1) (17.4.5)
Both equations (17.4.1) and (17.4.4) are invariant under the replacement
x → −x. Thus the functions S and y must also be invariant, except possibly for an
overall scale factor. (Since the equations are linear, a constant multiple of a solution
is also a solution.) Because the solutions will be normalized, the scale factor can
only be ±1. If n − m is odd, there are an odd number of zeros in the interval (−1, 1).
Thus we must choose the antisymmetric solution y(−x) = −y(x) which has a zero
at x = 0. Conversely, if n − m is even we must have the symmetric solution. Thus
ymn (−x) = (−1)n−m ymn (x) (17.4.6)

774 Chapter 17. Two Point Boundary Value Problems
and similarly for Smn .
The boundary conditions on (17.4.4) require that y be regular at x = ±1. In
other words, near the endpoints the solution takes the form
y = a0 + a1 (1 − x2 ) + a2 (1 − x2 )2 + . . . (17.4.7)
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Substituting this expansion in equation (17.4.4) and letting x → 1, we ﬁnd that
µ − c2
a1 = − a0 (17.4.8)
4(m + 1)
Equivalently,
µ − c2
y (1) = y(1) (17.4.9)
2(m + 1)
A similar equation holds at x = −1 with a minus sign on the right-hand side.
The irregular solution has a different relation between function and derivative at
the endpoints.
Instead of integrating the equation from −1 to 1, we can exploit the symmetry
(17.4.6) to integrate from 0 to 1. The boundary condition at x = 0 is
y(0) = 0, n − m odd
(17.4.10)
y (0) = 0, n − m even
A third boundary condition comes from the fact that any constant multiple
of a solution y is a solution. We can thus normalize the solution. We adopt the
m
normalization that the function Smn has the same limiting behavior as Pn at x = 1:
lim (1 − x2 )−m/2 Smn (x; c) = lim (1 − x2 )−m/2 Pn (x)
m
(17.4.11)
x→1 x→1
Various normalization conventions in the literature are tabulated by Flammer [1].
Imposing three boundary conditions for the second-order equation (17.4.4)
turns it into an eigenvalue problem for λ or equivalently for µ. We write it in the
standard form by setting
y1 = y (17.4.12)
y2 = y (17.4.13)
y3 = µ (17.4.14)
Then
y1 = y2 (17.4.15)
1
y2 = 2x(m + 1)y2 − (y3 − c2 x2 )y1 (17.4.16)
1 − x2
y3 = 0 (17.4.17)

17.5 Automated Allocation of Mesh Points 783
17.5 Automated Allocation of Mesh Points
In relaxation problems, you have to choose values for the independent variable at the
mesh points. This is called allocating the grid or mesh. The usual procedure is to pick
a plausible set of values and, if it works, to be content. If it doesn’t work, increasing the
number of points usually cures the problem.
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
If we know ahead of time where our solutions will be rapidly varying, we can put more
grid points there and less elsewhere. Alternatively, we can solve the problem ﬁrst on a uniform
mesh and then examine the solution to see where we should add more points. We then repeat
the solution with the improved grid. The object of the exercise is to allocate points in such
a way as to represent the solution accurately.
It is also possible to automate the allocation of mesh points, so that it is done
“dynamically” during the relaxation process. This powerful technique not only improves
the accuracy of the relaxation method, but also (as we will see in the next section) allows
internal singularities to be handled in quite a neat way. Here we learn how to accomplish
the automatic allocation.
We want to focus attention on the independent variable x, and consider two alternative
reparametrizations of it. The ﬁrst, we term q; this is just the coordinate corresponding to the
mesh points themselves, so that q = 1 at k = 1, q = 2 at k = 2, and so on. Between any two
mesh points we have ∆q = 1. In the change of independent variable in the ODEs from x to q,
dy
=g (17.5.1)
dx
becomes
dy dx
=g (17.5.2)
dq dq
In terms of q, equation (17.5.2) as an FDE might be written
dx dx
yk − yk−1 − 1
2 g + g =0 (17.5.3)
dq dq
k k−1
or some related version. Note that dx/dq should accompany g. The transformation between
x and q depends only on the Jacobian dx/dq. Its reciprocal dq/dx is proportional to the
density of mesh points.
Now, given the function y(x), or its approximation at the current stage of relaxation,
we are supposed to have some idea of how we want to specify the density of mesh points.
For example, we might want dq/dx to be larger where y is changing rapidly, or near to the
boundaries, or both. In fact, we can probably make up a formula for what we would like
dq/dx to be proportional to. The problem is that we do not know the proportionality constant.
That is, the formula that we might invent would not have the correct integral over the whole
range of x so as to make q vary from 1 to M , according to its deﬁnition. To solve this problem
we introduce a second reparametrization Q(q), where Q is a new independent variable. The
relation between Q and q is taken to be linear, so that a mesh spacing formula for dQ/dx
differs only in its unknown proportionality constant. A linear relation implies
d2 Q
=0 (17.5.4)
dq2
or, expressed in the usual manner as coupled ﬁrst-order equations,
dQ(x) dψ
=ψ =0 (17.5.5)
dq dq
where ψ is a new intermediate variable. We add these two equations to the set of ODEs
being solved.
Completing the prescription, we add a third ODE that is just our desired mesh-density
function, namely
dQ dQ dq
φ(x) = = (17.5.6)
dx dq dx