772 Chapter 17.Two Point Boundary Value Problems for (j-jz1;j<-jz2;j++){ Loop over columns to be zeroed. for(1=jm1:1<=jm2:1+) Loop over columns altered. vx-c[ic][l+loff][kc]; for(i=1z1;i<=1z2;i+)s[1][1]-=s[1][j]*vx; Loop over rows. vx=c[ic][jcf][kc]; for(i=iz1;1<=iz2;1++)s[i][jmf]-=s[1][j]*vx; Plus final element. 1c+=1; "Algebraically Difficult"Sets of Differential Equations 8= Relaxation methods allow you to take advantage of an additional opportunity that,while not obvious,can speed up some calculations enormously.It is not necessary that the set of variables y.correspond exactly with the dependent variables of the original differential equations.They can be related to those variables through algebraic equations.Obviously,it is necessary only that the solution variables allow us to evaluate the functions y,g,B,C that are used to construct the FDEs from the ODEs.In some problems g depends on functions of y that are known only implicitly,so that iterative solutions are necessary to evaluate functions in the ODEs.Often one can dispense with this "internal"nonlinear problem by defining a new set of variables from which both y,g and the boundary conditions can be obtained directly.A typical example occurs in physical problems where the equations require solution of a complex equation of state that can be expressed in more convenient terms using variables s是%ad的 2 Press. other than the original dependent variables in the ODE.While this approach is analogous to performing an analytic change of variables directly on the original ODEs,such an analytic transformation might be prohibitively complicated.The change of variables in the relaxation method is easy and requires no analytic manipulations. SCIENTIFIC CITED REFERENCES AND FURTHER READING: Eggleton,P.P.1971,Monthly Notices of the Royal Astronomical Society,vol.151,pp.351-364.[1] 6 Keller,H.B.1968,Numerical Methods for Two-Point Boundary-Value Problems (Waltham,MA: Blaisdell). Kippenhan,R.,Weigert,A.,and Hofmeister,E.1968,in Methods in Computational Physics, vol.7 (New York:Academic Press),pp.129ff. Numerica 10621 17.4 A Worked Example:Spheroidal Harmonics 431 Recipes The best way to understand the algorithms of the previous sections is to see (outside them employed to solve an actual problem.As a sample problem,we have selected the computation of spheroidal harmonics.(The more common name is spheroidal North angle functions,but we prefer the explicit reminder of the kinship with spherical harmonics.)We will show how to find spheroidal harmonics,first by the method of relaxation ($17.3),and then by the methods of shooting (817.1)and shooting to a fitting point ($17.2). Spheroidal harmonics typically arise when certain partial differential equations are solved by separation of variables in spheroidal coordinates.They satisfy the following differential equation on the interval-1 <x<1: a-+(-x-)s=0 (17.4.1)
772 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). for (j=jz1;j<=jz2;j++) { Loop over columns to be zeroed. for (l=jm1;l<=jm2;l++) { Loop over columns altered. vx=c[ic][l+loff][kc]; for (i=iz1;i<=iz2;i++) s[i][l] -= s[i][j]*vx; Loop over rows. } vx=c[ic][jcf][kc]; for (i=iz1;i<=iz2;i++) s[i][jmf] -= s[i][j]*vx; Plus final element. ic += 1; } } “Algebraically Difficult” Sets of Differential Equations Relaxation methods allow you to take advantage of an additional opportunity that, while not obvious, can speed up some calculations enormously. It is not necessary that the set of variables yj,k correspond exactly with the dependent variables of the original differential equations. They can be related to those variables through algebraic equations. Obviously, it is necessary only that the solution variables allow us to evaluate the functions y, g,B, C that are used to construct the FDEs from the ODEs. In some problems g depends on functions of y that are known only implicitly, so that iterative solutions are necessary to evaluate functions in the ODEs. Often one can dispense with this “internal” nonlinear problem by defining a new set of variables from which both y, g and the boundary conditions can be obtained directly. A typical example occurs in physical problems where the equations require solution of a complex equation of state that can be expressed in more convenient terms using variables other than the original dependent variables in the ODE. While this approach is analogous to performing an analytic change of variables directly on the original ODEs, such an analytic transformation might be prohibitively complicated. The change of variables in the relaxation method is easy and requires no analytic manipulations. CITED REFERENCES AND FURTHER READING: Eggleton, P.P. 1971, Monthly Notices of the Royal Astronomical Society, vol. 151, pp. 351–364. [1] Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell). Kippenhan, R., Weigert, A., and Hofmeister, E. 1968, in Methods in Computational Physics, vol. 7 (New York: Academic Press), pp. 129ff. 17.4 A Worked Example: Spheroidal Harmonics The best way to understand the algorithms of the previous sections is to see them employed to solve an actual problem. As a sample problem, we have selected the computation of spheroidal harmonics. (The more common name is spheroidal angle functions, but we prefer the explicit reminder of the kinship with spherical harmonics.) We will show how to find spheroidal harmonics, first by the method of relaxation (§17.3), and then by the methods of shooting (§17.1) and shooting to a fitting point (§17.2). Spheroidal harmonics typically arise when certain partial differential equations are solved by separation of variables in spheroidal coordinates. They satisfy the following differential equation on the interval −1 ≤ x ≤ 1: d dx (1 − x2) dS dx + λ − c2x2 − m2 1 − x2 S =0 (17.4.1)
17.4 A Worked Example:Spheroidal Harmonics 773 Here m is an integer,c is the"oblateness parameter,"and A is the eigenvalue.Despite the notation,c2 can be positive or negative.For c2>0 the functions are called “prolate,.”while if c20 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/2y (17.4.3) We then find from(17.4.1)that y satisfies the equation Numerica 10621 (1-x2) 2-2m+1)z Py dx +(4-c2x2)y=0 (17.4.4) where μ≡入-m(m+1) (17.4.5) Both equations (17.4.1)and (17.4.4)are invariant under the replacement --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.Ifn-m is odd,there are an odd number ofzeros in the interval (-1,1). Thus we must choose the antisymmetric solution y(-)=-y(x)which has a zero at =0.Conversely,if n-m is even we must have the symmetric solution.Thus ymn(-x)=(-1)-mymn(x) (17.4.6
17.4 A Worked Example: Spheroidal Harmonics 773 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 first the spherical case, where c = 0, we recognize the differential equation for Legendre functions P m n (x). In this case the eigenvalues are λmn = n(n + 1), n = m, m + 1,... . The integer n labels successive eigenvalues for fixed 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 fixed m, n = m, m + 1,... labels the successive eigenvalues. The computation of λmn(c) and Smn(x; c) traditionally has been quite difficult. 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 first 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) α∞ k=0 ak(1 ± x) k (17.4.2) in equation (17.4.1), we find 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/2y (17.4.3) We then find from (17.4.1) that y satisfies the equation (1 − x2) d2y dx2 − 2(m + 1)x dy dx + (µ − c2x2)y =0 (17.4.4) 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−mymn(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) Substituting this expansion in equation(17.4.4)and letting x-1,we find that 4-c2 a1=-4m+0 (17.4.8) Equivalently, 离 0-0 (17.4.9) A similar equation holds at z=-1 with a minus sign on the right-hand side. RECIPES I 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 Press. (17.4.6)to integrate from 0 to 1.The boundary condition at x=0 is R y(0)=0,n-m odd Progra (0)=0,n-m even (17.4.10) OF SCIENTIFIC( 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 normalization that the function Smn has the same limiting behavior as Pm at =1: lim(1-22)-m/2Smn(;c)=lim(1-22)-m/2p() (17.4.11) 10621 Various normalization conventions in the literature are tabulated by Flammer [1]. Numerical Imposing three boundary conditions for the second-order equation(17.4.4) 431 turns it into an eigenvalue problem for A or equivalently for u.We write it in the Recipes standard form by setting 1=y (17.4.12) 2=y (17.4.13) 3=4 (17.4.14) Then =2 (17.4.15) 1 =12 2=(m+1)v2-(u3-)] (17.4.16) 5=0 (17.4.17)
774 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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) Substituting this expansion in equation (17.4.4) and letting x → 1, we find that a1 = − µ − c2 4(m + 1)a0 (17.4.8) Equivalently, y (1) = µ − c2 2(m + 1)y(1) (17.4.9) 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 y (0) = 0, n − m even (17.4.10) 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 normalization that the function Smn has the same limiting behavior as P m n at x = 1: limx→1 (1 − x2) −m/2Smn(x; c) = limx→1 (1 − x2) −m/2P m n (x) (17.4.11) 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 y 1 = y2 (17.4.15) y 2 = 1 1 − x2 2x(m + 1)y2 − (y3 − c2x2)y1 (17.4.16) y 3 =0 (17.4.17)
17.4 A Worked Example:Spheroidal Harmonics 775 The boundary condition at z=0 in this notation is y1=0,n-m odd (17.4.18) y2 =0,n-m even At x =1 we have two conditions: 2 3-c2 2(m+1) (17.4.19) h-1im(1-x2)m/2P(a)= (-1)m(m+m)! 2mm!(n-m)! 三Y (17.4.20) We are now ready to illustrate the use of the methods of previous sections on this problem. Relaxation g经 If we just want a few isolated values of A or S,shooting is probably the quickest University 2 method.However,if we want values for a large sequence of values of c,relaxation is better.Relaxation rewards a good initial guess with rapid convergence,and the Press. THE previous solution should be a good initial guess if c is changed only slightly. For simplicity,we choose a uniform grid on the interval 0 <z<1.For a total of M mesh points,we have Program 1 h=M-1 (17.4.21) OF SCIENTIFIC( xk=(k-1)h, k=1,2,.,M (17.4.22) At interior points k =2,3,...,M,equation (17.4.15)gives COMPUTING (ISBN h 19200 E1.k =y1.k-y1.k-1-v2k+y2k-1) (17.4.23) 10-521 Equation (17.4.16)gives 43106 E2.k=2,k-2,k-1-Bk [k+k-m+12k+2k-1-ank+1k- (17.4.24) (outside + North Software. where a4=%k+3-1-2(%+xk-1)2 (17.4.25) 2 4 h 脉=1-+k- (17.4.26) Finally,equation (17.4.17)gives E3.k=3,k-3,k-1 (17.4.27)
17.4 A Worked Example: Spheroidal Harmonics 775 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The boundary condition at x = 0 in this notation is y1 = 0, n − m odd y2 = 0, n − m even (17.4.18) At x = 1 we have two conditions: y2 = y3 − c2 2(m + 1)y1 (17.4.19) y1 = limx→1 (1 − x2) −m/2P m n (x) = (−1)m(n + m)! 2mm!(n − m)! ≡ γ (17.4.20) We are now ready to illustrate the use of the methods of previous sections on this problem. Relaxation If we just want a few isolated values of λ or S, shooting is probably the quickest method. However, if we want values for a large sequence of values of c, relaxation is better. Relaxation rewards a good initial guess with rapid convergence, and the previous solution should be a good initial guess if c is changed only slightly. For simplicity, we choose a uniform grid on the interval 0 ≤ x ≤ 1. For a total of M mesh points, we have h = 1 M − 1 (17.4.21) xk = (k − 1)h, k = 1, 2,...,M (17.4.22) At interior points k = 2, 3,...,M, equation (17.4.15) gives E1,k = y1,k − y1,k−1 − h 2 (y2,k + y2,k−1) (17.4.23) Equation (17.4.16) gives E2,k = y2,k − y2,k−1 − βk × (xk + xk−1)(m + 1)(y2,k + y2,k−1) 2 − αk (y1,k + y1,k−1) 2 (17.4.24) where αk = y3,k + y3,k−1 2 − c2(xk + xk−1)2 4 (17.4.25) βk = h 1 − 1 4 (xk + xk−1)2 (17.4.26) Finally, equation (17.4.17) gives E3,k = y3,k − y3,k−1 (17.4.27)
776 Chapter 17.Two Point Boundary Value Problems Now recall that the matrix of partial derivatives Si of equation(17.3.8)is defined so that i labels the equation and j the variable.In our case,j runs from I to 3 for yi at k-1 and from 4 to 6 for yi at k.Thus equation(17.4.23)gives h S1,1=-1, 512=-21 S1.3=0 (17.4.28) S1,4=1, 515=-2 S1,6=0 Similarly equation (17.4.24)yields S2,1=ak3k/2, S2,2=-1-k(xk+xk-1)(m+1)/2, S2,3=Bk(1,k+1,k-1)/4 S2,4=S2.1, S2,5=2+S2,2 S2,6=S2.3 (17.4.29) while from equation (17.4.27)we find RECIPES S3.1=0, S32=0, S3,3=-1 2 (17.4.30) S3.4=0, S3.5=0, S3.6=1 At z =0 we have the boundary condition Press. E3.1= ∫h,1,n-modd (17.4.31) Programs 22.1,n-m even Recall the convention adopted in the solvde routine that for one boundary condition at k =1 only S3.i can be nonzero.Also,j takes on the values 4 to 6 since the 6 boundary condition involves only yk,not yk-1.Accordingly,the only nonzero values of S3.j at x =0 are S3.4=1, n-m odd (17.4.32) S3.5=1, n-m even At z =1 we have 然 Numerical 10621 431 Recipes E1,M+1=欢.M- 3M-c2 2m+万h,M (17.4.33) (outside E2,M+1=1,M-Y (17.4.34) Thus 3M-c2 51,4=-2m+丁可 S1,5=1, 1,M S1,6=-2m+可 (17.4.35) S2,4=1, S2,5=0, S2.6=0 (17.4.36 Here now is the sample program that implements the above algorithm.We need a main program,sfroid,that calls the routine solvde,and we must supply the function difeq called by solvde.For simplicity we choose an equally spaced
776 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Now recall that the matrix of partial derivatives Si,j of equation (17.3.8) is defined so that i labels the equation and j the variable. In our case, j runs from 1 to 3 for yj at k − 1 and from 4 to 6 for yj at k. Thus equation (17.4.23) gives S1,1 = −1, S1,2 = −h 2 , S1,3 = 0 S1,4 = 1, S1,5 = −h 2 , S1,6 = 0 (17.4.28) Similarly equation (17.4.24) yields S2,1 = αkβk/2, S2,2 = −1 − βk(xk + xk−1)(m + 1)/2, S2,3 = βk(y1,k + y1,k−1)/4 S2,4 = S2,1, S2,5 =2+ S2,2, S2,6 = S2,3 (17.4.29) while from equation (17.4.27) we find S3,1 = 0, S3,2 = 0, S3,3 = −1 S3,4 = 0, S3,5 = 0, S3,6 = 1 (17.4.30) At x = 0 we have the boundary condition E3,1 = y1,1, n − m odd y2,1, n − m even (17.4.31) Recall the convention adopted in the solvde routine that for one boundary condition at k = 1 only S3,j can be nonzero. Also, j takes on the values 4 to 6 since the boundary condition involves only y k, not yk−1. Accordingly, the only nonzero values of S3,j at x = 0 are S3,4 = 1, n − m odd S3,5 = 1, n − m even (17.4.32) At x = 1 we have E1,M+1 = y2,M − y3,M − c2 2(m + 1) y1,M (17.4.33) E2,M+1 = y1,M − γ (17.4.34) Thus S1,4 = −y3,M − c2 2(m + 1) , S1,5 = 1, S1,6 = − y1,M 2(m + 1) (17.4.35) S2,4 = 1, S2,5 = 0, S2,6 = 0 (17.4.36) Here now is the sample program that implements the above algorithm. We need a main program, sfroid, that calls the routine solvde, and we must supply the function difeq called by solvde. For simplicity we choose an equally spaced
17.4 A Worked Example:Spheroidal Harmonics 777 mesh of m=41 points,that is,h=.025.As we shall see,this gives good accuracy for the eigenvalues up to moderate values of n-m. Since the boundary condition at z =0 does not involve y if n-m is even, we have to use the indexv feature of solvde.Recall that the value of indexv [j] describes which column of s[i][j]the variable y[j]has been put in.If n-m is even,we need to interchange the columns for y and y2 so that there is not a zero pivot element in s[i][j] The program prompts for values of m and n.It then computes an initial guess for y based on the Legendre function Pm.It next prompts for c2,solves for y, prompts for c2,solves for y using the previous values as an initial guess,and so on. #include #include nted for 18881992 #include "nrutil.h" #define NE 3 100 (including this one) #define M 41 #define NB 1 872 #define NSI NE 7423 users from NUMERICAL RECIPES IN #define NYJ NE #define NYK M #define NCI NE #define NCJ (NE-NB+1) #define NCK (M+1) (North America to to sto make one paper Cambridge University Press. THE #define NSJ (2*NE+1) 是 ART int mm,n,mpt=M; float h,c2=0.0,anorm,x[M+1]; Global variables communicating with difeq. int main(void)/*Program sfroid * strictly prohibited Programs Sample program using solvde.Computes eigenvalues of spheroidal harmonics Smn(z:c)for m 20 and n m.In the program,m is mm,c2 is c2,and y of equation (17.4.20)is anorm. to dir Copyright (C) float plgndr(int 1,int m,float x); void solvde(int itmax,float conv,float slowc,float scalv[], int indexv[],int ne,int nb,int m,float **y,float ***c, float **s); int i,itmax,k,indexv [NE+1]; OF SCIENTIFIC COMPUTING (ISBN 0-521 float conv,deriv,fac1,fac2,q1,slowc,scalv[NE+1]; f10at米*y,*s,米*C; y=matrix(1,NYJ,1,NYK): s=matrix(1,NSI,1,NSJ); 1988-1992 by Numerical Recipes c=f3tensor(1,NCI,1,NCJ,1,NCK); -43108 itmax=100; conv=5.0e-6; s1owc=1.0; h=1.0/(M-1); Software. printf("\nenter m n\n"); scanf ("%d %d",&mm,&n); if (n+mm &1){ No interchanges necessary (outside North America) indexv [1]=1; indexv [2]=2; indexv [3]=3; else Interchange y1 and y2. indexv[1]=2; indexv[2]=1; indexv [3]=3; anorm=1.0; Compute Y. if (mm){ qi=n;
17.4 A Worked Example: Spheroidal Harmonics 777 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). mesh of m = 41 points, that is, h = .025. As we shall see, this gives good accuracy for the eigenvalues up to moderate values of n − m. Since the boundary condition at x = 0 does not involve y 1 if n − m is even, we have to use the indexv feature of solvde. Recall that the value of indexv[j] describes which column of s[i][j] the variable y[j] has been put in. If n − m is even, we need to interchange the columns for y 1 and y2 so that there is not a zero pivot element in s[i][j]. The program prompts for values of m and n. It then computes an initial guess for y based on the Legendre function P m n . It next prompts for c2, solves for y, prompts for c2, solves for y using the previous values as an initial guess, and so on. #include #include #include "nrutil.h" #define NE 3 #define M 41 #define NB 1 #define NSI NE #define NYJ NE #define NYK M #define NCI NE #define NCJ (NE-NB+1) #define NCK (M+1) #define NSJ (2*NE+1) int mm,n,mpt=M; float h,c2=0.0,anorm,x[M+1]; Global variables communicating with difeq. int main(void) /* Program sfroid */ Sample program using solvde. Computes eigenvalues of spheroidal harmonics Smn(x; c) for m ≥ 0 and n ≥ m. In the program, m is mm, c2 is c2, and γ of equation (17.4.20) is anorm. { float plgndr(int l, int m, float x); void solvde(int itmax, float conv, float slowc, float scalv[], int indexv[], int ne, int nb, int m, float **y, float ***c, float **s); int i,itmax,k,indexv[NE+1]; float conv,deriv,fac1,fac2,q1,slowc,scalv[NE+1]; float **y,**s,***c; y=matrix(1,NYJ,1,NYK); s=matrix(1,NSI,1,NSJ); c=f3tensor(1,NCI,1,NCJ,1,NCK); itmax=100; conv=5.0e-6; slowc=1.0; h=1.0/(M-1); printf("\nenter m n\n"); scanf("%d %d",&mm,&n); if (n+mm & 1) { No interchanges necessary. indexv[1]=1; indexv[2]=2; indexv[3]=3; } else { Interchange y1 and y2. indexv[1]=2; indexv[2]=1; indexv[3]=3; } anorm=1.0; Compute γ. if (mm) { q1=n;
778 Chapter 17.Two Point Boundary Value Problems for (i=1;iscalv[1]?y[2]0:sca1v[1]); 1.800 sca1v[3]=(y[3]0>1.0?y[3]0:1.0); for(;;){ 872 printf("\nEnter c**2 or 999 to end.In"); to any Cambridge from NUMERICAL RECIPESI scanf("%f",&c2); 1f(c2==999)[ free_f3tensor(c,1,NCI,1,NCJ,1,NCK); free_matrix(s,1,NSI,1,NSJ); free_matrix(y,1,NYJ,1,NYK); THE return 0; server computer, (North America one paper University Press. ART solvde(itmax,conv,slowc,scalv,indexv,NE,NB,M,y,c,s); printf("\n %s %2d %s %2d %s %7.3f %s %10.6f\n", m=",mm,"n=",n,"c*2=",c2, "1amda=",y[3][1]+mm*(mm+1)); strictly proh Programs Return for another value of c2. extern int mm,n,mpt; Defined in sfroid. 188 extern float h,c2,anorm,x[]; OF SCIENTIFIC COMPUTING(ISBN 1920 void difeg(int k,int ki,int k2,int isf,int is1,int isf,int indexv[], int ne,f1oat米*s,f1oat来*y) Returns matrix s for solvde. float temp,tempi,temp2; Numerical Recipes 021 43108 1f(k==k1){ Boundary condition at first point. if (n+mm 1){ g[3][3+indexv[1]]=1.0; Equation (17.4.32) (outside 膜 s[3][3+indexv[2]]=0.0; s[3][3+indexv[3]]=0.0; oftware. s[3][jsf]=y[1][1]; Equation (17.4.31). else s[3][3+indexv[1]]=0.0; Equation (17.4.32) s[3][3+indexv[2]]=1.0; s[3][3+indexv[3]]=0.0; s[3][jsf]=y[2][1]; Equation (17.4.31) else if (k>k2){ Boundary conditions at last point. s[1][3+indexv[1]]=-(y[3][mpt]-c2)/(2.0*(mm+1.0)); (17.4.35) s[1][3+indexv[2]=1.0; s[1][3+indexv[3]]=-y[1][mpt]/(2.0*(mm+1.0)); s[1][jsf]=y[2][mpt]-(y[3][mpt]-c2)*y[1][mpt]/(2.0*(mm+1.0));(17.4.33). s[2][3+indexv[1]]=1.0; Equation (17.4.36)
778 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). for (i=1;i scalv[1] ? y[2][M] : scalv[1]); scalv[3]=(y[3][M] > 1.0 ? y[3][M] : 1.0); for (;;) { printf("\nEnter c**2 or 999 to end.\n"); scanf("%f",&c2); if (c2 == 999) { free_f3tensor(c,1,NCI,1,NCJ,1,NCK); free_matrix(s,1,NSI,1,NSJ); free_matrix(y,1,NYJ,1,NYK); return 0; } solvde(itmax,conv,slowc,scalv,indexv,NE,NB,M,y,c,s); printf("\n %s %2d %s %2d %s %7.3f %s %10.6f\n", "m =",mm," n =",n," c**2 =",c2, " lamda =",y[3][1]+mm*(mm+1)); } Return for another value of c2. } extern int mm,n,mpt; Defined in sfroid. extern float h,c2,anorm,x[]; void difeq(int k, int k1, int k2, int jsf, int is1, int isf, int indexv[], int ne, float **s, float **y) Returns matrix s for solvde. { float temp,temp1,temp2; if (k == k1) { Boundary condition at first point. if (n+mm & 1) { s[3][3+indexv[1]]=1.0; Equation (17.4.32). s[3][3+indexv[2]]=0.0; s[3][3+indexv[3]]=0.0; s[3][jsf]=y[1][1]; Equation (17.4.31). } else { s[3][3+indexv[1]]=0.0; Equation (17.4.32). s[3][3+indexv[2]]=1.0; s[3][3+indexv[3]]=0.0; s[3][jsf]=y[2][1]; Equation (17.4.31). } } else if (k > k2) { Boundary conditions at last point. s[1][3+indexv[1]] = -(y[3][mpt]-c2)/(2.0*(mm+1.0)); (17.4.35). s[1][3+indexv[2]]=1.0; s[1][3+indexv[3]] = -y[1][mpt]/(2.0*(mm+1.0)); s[1][jsf]=y[2][mpt]-(y[3][mpt]-c2)*y[1][mpt]/(2.0*(mm+1.0)); (17.4.33). s[2][3+indexv[1]]=1.0; Equation (17.4.36)
17.4 A Worked Example:Spheroidal Harmonics 779 s[2][3+indexv[2]]=0.0: s[2][3+indexv[3]]=0.0; s[2][jsf]=y[1][mpt]-anorm; Equation (17.4.34) else Interior point. s[1][indexv[1]]=-1.0: Equation (17.4.28). s[1][indexv[2]]=-0.5*h; s[1][indexv[3]=0.0; s[1][3+indexv[1]]=1.0: s[1][3+indexv[2]]=-0.5wh; s[1][3+indexv[3]]=0.0; temp1=x[k]+x[k-1]; temp=h/(1.0-templ*temp1*0.25): http://www.n read able files temp2=0.5*(y[3][k]+y[3][k-1])-c2*0.25*temp1*temp1; s[2][indexv [1]]=temp*temp2*0.5; Equation (17.4.29) s[2][indexv[2]=-1.0-0.5*temp*(mm+1.0)*temp1; 83g s[2][indexv[3]]=0.25*temp*(y[1][x]+y[1][k-1]); (including this one) granted for 19881992 s[2][3+indexv [1]]=s[2][indexv[1]]; s[2][3+indexv[2]]=2.0+s[2][indexv[2]]; s[2][3+indexv [3]]=s [2][indexv [3]]; s[3][indexv[1]]=0.0; Equation (17.4.30). 111800-672 s[3][indexv[2]]=0.0; to any /Cambridge s[3][indexv[3]=-1.0; from NUMERICAL RECIPES IN s[3][3+indexv[1]]=0.0: s[3][3+indexv[2]]=0.0; s[3][3+1 ndexv[3]]=1.0; (Nort s[1][jsf]=y[1][k]-y[1][k-1]-0.5*h*(y[2][k]+y[2][k-1]): (17.4.23)」 s[2][jsf]=y[2][k]-y[2]k-1]-temp*((xk]+x[k-1]) (17.4.24). America server computer, to make one paper UnN电.t THE *0.5*(mm+1.0)*(y[2][k]+y[2][k-1])-temp2 ART *0.5*(y[1]]+y[1][k-1]); s[3][jsf]=y[3][k]-y[3][k-1]; Equation (17.4.27). Programs 2 You can run the program and check it against values ofAmn(c)given in the tables at the back of Flammer's book [1]or in Table 21.1 of Abramowitz and Stegun [21. to dir Typically it converges in about 3 iterations.The table below gives a few comparisons. Selected Output of sfroid m n 3 入exact 入sfroid 2 2 0.1 6.01427 6.01427 1.0 6.14095 6.14095 之究 Further reproduction,or 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 10-:6211 43108 4.0 6.54250 6.54253 1.0 30.4361 30.4372 16.0 36.9963 37.0135 Software. 11 -1.0 131.560 131.554 (outside North America). Shooting To solve the same problem via shooting(817.1),we supply a function derivs that implements equations(17.4.15)(17.4.17).We will integrate the equations over the range-1 <z<0.We provide the function load which sets the eigenvalue y3 to its current best estimate,v[].It also sets the boundary values of y and y2 using equations(17.4.20)and(17.4.19)(with a minus sign corresponding to x=-1).Note that the boundary condition is actually applied a distance dx from
17.4 A Worked Example: Spheroidal Harmonics 779 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). s[2][3+indexv[2]]=0.0; s[2][3+indexv[3]]=0.0; s[2][jsf]=y[1][mpt]-anorm; Equation (17.4.34). } else { Interior point. s[1][indexv[1]] = -1.0; Equation (17.4.28). s[1][indexv[2]] = -0.5*h; s[1][indexv[3]]=0.0; s[1][3+indexv[1]]=1.0; s[1][3+indexv[2]] = -0.5*h; s[1][3+indexv[3]]=0.0; temp1=x[k]+x[k-1]; temp=h/(1.0-temp1*temp1*0.25); temp2=0.5*(y[3][k]+y[3][k-1])-c2*0.25*temp1*temp1; s[2][indexv[1]]=temp*temp2*0.5; Equation (17.4.29). s[2][indexv[2]] = -1.0-0.5*temp*(mm+1.0)*temp1; s[2][indexv[3]]=0.25*temp*(y[1][k]+y[1][k-1]); s[2][3+indexv[1]]=s[2][indexv[1]]; s[2][3+indexv[2]]=2.0+s[2][indexv[2]]; s[2][3+indexv[3]]=s[2][indexv[3]]; s[3][indexv[1]]=0.0; Equation (17.4.30). s[3][indexv[2]]=0.0; s[3][indexv[3]] = -1.0; s[3][3+indexv[1]]=0.0; s[3][3+indexv[2]]=0.0; s[3][3+indexv[3]]=1.0; s[1][jsf]=y[1][k]-y[1][k-1]-0.5*h*(y[2][k]+y[2][k-1]); (17.4.23). s[2][jsf]=y[2][k]-y[2][k-1]-temp*((x[k]+x[k-1]) (17.4.24). *0.5*(mm+1.0)*(y[2][k]+y[2][k-1])-temp2 *0.5*(y[1][k]+y[1][k-1])); s[3][jsf]=y[3][k]-y[3][k-1]; Equation (17.4.27). } } You can run the program and check it against values of λ mn(c) given in the tables at the back of Flammer’s book [1] or in Table 21.1 of Abramowitz and Stegun [2]. Typically it converges in about 3 iterations. The table below gives a few comparisons. Selected Output of sfroid mn c2 λexact λsfroid 22 0.1 6.01427 6.01427 1.0 6.14095 6.14095 4.0 6.54250 6.54253 25 1.0 30.4361 30.4372 16.0 36.9963 37.0135 4 11 −1.0 131.560 131.554 Shooting To solve the same problem via shooting (§17.1), we supply a function derivs that implements equations (17.4.15)–(17.4.17). We will integrate the equations over the range −1 ≤ x ≤ 0. We provide the function load which sets the eigenvalue y3 to its current best estimate, v[1]. It also sets the boundary values of y 1 and y2 using equations (17.4.20) and (17.4.19) (with a minus sign corresponding to x = −1). Note that the boundary condition is actually applied a distance dx from
780 Chapter 17.Two Point Boundary Value Problems the boundary to avoid having to evaluate y2 right on the boundary.The function score follows from equation (17.4.18). #include #include "nrutil.h' #define N2 1 int m,n; Communicates with load,score,and derivs. float c2,dx,gmmai int nvar; Communicates with shoot. float x1,x2; 83g int main(void)Program sphoot * Sample program using shoot.Computes eigenvalues of spheroidal harmonics Smn(z;c)for granted for 19881992 m >0 and n>m.Note how the routine vecfunc for newt is provided by shoot ($17.1). 100 void newt(float x[],int n,int *check, void(*vecfunc)(int,float[☐,f1oatd)): void shoot(int n,float v[],float f[); to any from NUMERICAL RECIPES IN int check,i; float q1,*v; v=vector(1,N2); (North America 州bMe se to make one paper by Cambridge University Press. dx=1.0e-4; Avoid evaluating derivatives exactly at x=-1. THE nvar=3; Number of equations. for(;;)[ printf("input m,n,c-squared\n"); 是 ART if (scanf("%d %d f",&m,&n,&c2)=EOF)break; if (n m II m 0)continue; gmma=1.0; Compute y of equation (17.4.20) q1=n; st st Programs for(1=1;i<=m;1++)gmma*=-0.5*(n+i)*(q1--/1); email v[1]=n*(n+1)-m*(m+1)+c2/2.0; Initial guess for eigenvalue Copyright (C) x1=-1.0+dx; Set range of integration. x2=0.0: newt (v,N2,&check,shoot); Find v that zeros function f in score. if (check){ printf("shoot failed;bad initial guess\n"); else OF SCIENTIFIC COMPUTING(ISBN 0-521- printf("\tmu(m,n)\n"); printf("%12.6f\n",v[1]) @cambridge.org 1988-1992 by Numerical Recipes free_vector(v,1,N2); -431085 return 0; void load(float x1,float v[],float y]) Software. Supplies starting values for integration at x =-1 +dr. float y1 (n-m&1 -gmmagmma); (outside North America) y[3]=v[1]; y[2]=-(y[3]-c2)*y1/(2*(m+1)); y[1]=y1+y[2]*dx; void score(float xf,float y[],float f[) Tests whether boundary condition at x =0 is satisfied. f[1]=(n-m&1?y[1]:y[2]);
780 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). the boundary to avoid having to evaluate y 2 right on the boundary. The function score follows from equation (17.4.18). #include #include "nrutil.h" #define N2 1 int m,n; Communicates with load, score, and derivs. float c2,dx,gmma; int nvar; Communicates with shoot. float x1,x2; int main(void) /* Program sphoot */ Sample program using shoot. Computes eigenvalues of spheroidal harmonics Smn(x; c) for m ≥ 0 and n ≥ m. Note how the routine vecfunc for newt is provided by shoot (§17.1). { void newt(float x[], int n, int *check, void (*vecfunc)(int, float [], float [])); void shoot(int n, float v[], float f[]); int check,i; float q1,*v; v=vector(1,N2); dx=1.0e-4; Avoid evaluating derivatives exactly at x = −1. nvar=3; Number of equations. for (;;) { printf("input m,n,c-squared\n"); if (scanf("%d %d %f",&m,&n,&c2) == EOF) break; if (n < m || m < 0) continue; gmma=1.0; Compute γ of equation (17.4.20). q1=n; for (i=1;i<=m;i++) gmma *= -0.5*(n+i)*(q1--/i); v[1]=n*(n+1)-m*(m+1)+c2/2.0; Initial guess for eigenvalue. x1 = -1.0+dx; Set range of integration. x2=0.0; newt(v,N2,&check,shoot); Find v that zeros function f in score. if (check) { printf("shoot failed; bad initial guess\n"); } else { printf("\tmu(m,n)\n"); printf("%12.6f\n",v[1]); } } free_vector(v,1,N2); return 0; } void load(float x1, float v[], float y[]) Supplies starting values for integration at x = −1 + dx. { float y1 = (n-m&1? -gmma : gmma); y[3]=v[1]; y[2] = -(y[3]-c2)*y1/(2*(m+1)); y[1]=y1+y[2]*dx; } void score(float xf, float y[], float f[]) Tests whether boundary condition at x = 0 is satisfied. { f[1]=(n-m & 1 ? y[1] : y[2]); }
17.4 A Worked Example:Spheroidal Harmonics 781 void derivs(float x,float y[],float dydx[]) Evaluates derivatives for odeint. dydx[1]=y[2]: dydx[2]=(2.0*x*(m+1.0)*y[2]-(y[3]-c2*x*x)*y[1])/(1.0-x*x); dydx[3]=0.0; Shooting to a Fitting Point For variety we illustrate shootf from $17.2 by integrating over the whole range 8 -1+dz #include ART #include "nrutil.h" #define N1 2 #define N2 1 Programs #define NTOT (N1+N2) #define DXX 1.0e-4 int m,n; Communicates with load1,load2,score, float c2,dx,gmma; and derivs. int nn2,nvar; Communicates with shootf. float x1,x2,xf; 19881992 SCIENTIFIC COMPUTING(ISBN int main(void)/*Program sphfpt * Sample program using shootf.Computes eigenvalues of spheroidal harmonics Smn(z;c)for m >0 and n m.Note how the routine vecfunc for newt is provided by shootf ($17.2). Numerical Recipes 10-621 The routine derivs is the same as for sphoot. 43108 void newt(float x[],int n,int *check, void(*vecfunc)(int,f1oat☐,f1oat0)); void shootf(int n,float v[],float ) int check,i; (outside float q1,*v1,*v2,*v; Software. v=vector(1,NTOT); v1=V; Amer v2=v[N2]; nvar=NTOT; Number of equations. nn2=N2; dx=DXX; Avoid evaluating derivatives exactly at z for (;;) 士1. printf("input m,n,c-squared\n"); if (scanf("%d %d %f",&m,&n,&c2)==EOF)break; if (n m Il m 0)continue; gmma=1.0; Compute y of equation (17.4.20). q1=n; for(1=1;1<=m;i+)gmma*=-0.5*(n+1)*(q1--/1);
17.4 A Worked Example: Spheroidal Harmonics 781 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). void derivs(float x, float y[], float dydx[]) Evaluates derivatives for odeint. { dydx[1]=y[2]; dydx[2]=(2.0*x*(m+1.0)*y[2]-(y[3]-c2*x*x)*y[1])/(1.0-x*x); dydx[3]=0.0; } Shooting to a Fitting Point For variety we illustrate shootf from §17.2 by integrating over the whole range −1 + dx ≤ x ≤ 1 − dx, with the fitting point chosen to be at x = 0. The routine derivs is identical to the one for shoot. Now, however, there are two load routines. The routine load1 for x = −1 is essentially identical to load above. At x = 1, load2 sets the function value y1 and the eigenvalue y3 to their best current estimates, v2[1] and v2[2], respectively. If you quite sensibly make your initial guess of the eigenvalue the same in the two intervals, then v1[1] will stay equal to v2[2] during the iteration. The function score simply checks whether all three function values match at the fitting point. #include #include #include "nrutil.h" #define N1 2 #define N2 1 #define NTOT (N1+N2) #define DXX 1.0e-4 int m,n; Communicates with load1, load2, score, float c2,dx,gmma; and derivs. int nn2,nvar; Communicates with shootf. float x1,x2,xf; int main(void) /* Program sphfpt */ Sample program using shootf. Computes eigenvalues of spheroidal harmonics Smn(x; c) for m ≥ 0 and n ≥ m. Note how the routine vecfunc for newt is provided by shootf (§17.2). The routine derivs is the same as for sphoot. { void newt(float x[], int n, int *check, void (*vecfunc)(int, float [], float [])); void shootf(int n, float v[], float f[]); int check,i; float q1,*v1,*v2,*v; v=vector(1,NTOT); v1=v; v2 = &v[N2]; nvar=NTOT; Number of equations. nn2=N2; dx=DXX; Avoid evaluating derivatives exactly at x = for (;;) { ±1. printf("input m,n,c-squared\n"); if (scanf("%d %d %f",&m,&n,&c2) == EOF) break; if (n < m || m < 0) continue; gmma=1.0; Compute γ of equation (17.4.20). q1=n; for (i=1;i<=m;i++) gmma *= -0.5*(n+i)*(q1--/i);