f17-4

764 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).17.4 AWorkedExample: SpheroidalHarmonics The best way to understand the algorithms of the previous sections is to see them employedto solve an actual problem. As a sample problem, we have selectedthe 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/bracketleftbigg (1−x2)dS dx/bracketrightbigg +/parenleftbigg λ−c2x2−m2 1−x2/parenrightbigg S=0 ( 17.4.1 ) Heremisaninteger, cisthe“oblatenessparameter,”and λistheeigenvalue. Despite the notation, c2can be positive or negative. For c2>0the functions are called “prolate,” while if c2<0they are called “oblate.” The equationhas singular points atx=±1andistobesolvedsubjecttotheboundaryconditionsthatthesolutionbe regularat x=±1. Onlyforcertainvaluesof λ,theeigenvalues,willthisbepossible. Ifweconsiderfirstthesphericalcase,where c=0,werecognizethedifferential equation for Legendre functions Pm n(x). In this case the eigenvalues are λmn = n(n+1 ),n=m, m +1,.... The integer nlabels successive eigenvalues for fixed m: When n=mwe have the lowest eigenvalue, and the corresponding eigenfunctionhas no nodes in the interval −1<x< 1; when n=m+1we have the nexteigenvalue,andthe eigenfunctionhas onenodeinside (−1,1); andso on. A similar situation holdsforthe generalcase c2/negationslash=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. Thecomputationof λmn(c)andSmn(x;c)traditionallyhasbeenquitedifficult. 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)α∞/summationdisplay k=0ak(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≥0since m→−mis a symmetry of the equation.) We get an equationthat is numericallymore tractable if we factor out this behavior.Accordingly we set S=( 1−x 2)m/2y (17.4.3 ) We then find from (17.4.1) that ysatisfies the equation (1−x2)d2y dx2−2(m+1 )xdy dx+(µ−c2x2)y=0 ( 17.4.4 ) 17.4A WorkedExample: SpheroidalHarmonics 765Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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 Sandymust also be invariant,except possibly for an overallscale factor. (Sincetheequationsarelinear,aconstantmultipleofa solution is also a solution.) Because the solutions will be normalized, the scale factor canonlybe ±1.I fn−misodd,thereareanoddnumberofzerosintheinterval (−1,1). Thus we must choosethe antisymmetricsolution y(−x)=−y(x)which has a zero atx=0. Conversely,if n−mis evenwe musthavethe symmetricsolution. Thus y mn(−x)=(−1)n−mymn(x)( 17.4.6 ) and similarly for Smn. The boundary conditions on (17.4.4) require that ybe regular at x=±1.I n other words, near the endpoints the solution takes the form y=a0+a1(1−x2)+a2(1−x2)2+... (17.4.7 ) Substituting this expansionin equation(17.4.4)and letting x→1, we find that a1=−µ−c2 4(m+1 )a0 (17.4.8 ) Equivalently, y/prime(1) =µ−c2 2(m+1 )y(1) ( 17.4.9 ) A similar equation holds at x=−1with a minus sign on the right-hand side. The irregular solution has a different relation between function and derivative atthe endpoints. Instead of integratingthe equation from −1to 1, we can exploit the symmetry (17.4.6) to integrate from 0 to 1. The boundary condition at x=0is y(0) = 0 ,n−modd y /prime(0) = 0 ,n−meven(17.4.10 ) A third boundary condition comes from the fact that any constant multiple of a solution yis a solution. We can thus normalize the solution. We adopt the normalizationthatthefunction Smnhasthesamelimitingbehavioras Pm natx=1: lim x→1(1−x2)−m/2Smn(x;c) = lim x→1(1−x2)−m/2Pm n(x)( 17.4.11 ) Variousnormalizationconventionsin the literature are tabulatedby Flammer [1]. 766 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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/prime(17.4.13 ) y3=µ (17.4.14 ) Then y/prime 1=y2 (17.4.15 ) y/prime 2=1 1−x2/bracketleftbig 2x(m+1 )y2−(y3−c2x2)y1/bracketrightbig (17.4.16 ) y/prime 3=0 ( 17.4.17 ) The boundary condition at x=0in this notation is y1=0,n−modd y2=0,n−meven(17.4.18 ) Atx=1we have two conditions: y2=y3−c2 2(m+1 )y1 (17.4.19 ) y1= lim x→1(1−x2)−m/2Pm 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 Ifwejustwantafewisolatedvaluesof λorS,shootingisprobablythequickest 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 cis changed only slightly. For simplicity, we choose a uniform grid on the interval 0≤x≤1. For a total of Mmesh 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 ) 17.4A WorkedExample: SpheroidalHarmonics 767Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Equation (17.4.16) gives E2,k=y2,k−y2,k−1−βk ×/bracketleftbigg(xk+xk−1)(m+1 ) ( y2,k+y2,k−1) 2−αk(y1,k+y1,k−1) 2/bracketrightbigg(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 ) Now recall that the matrix of partial derivatives Si,jof equation (17.3.8) is definedso that ilabels theequationand jthe variable. In ourcase, jrunsfrom1 to 3 for yjatk−1and from 4 to 6 for yjatk. Thus equation (17.4.23)gives S1,1=−1,S 1,2=−h 2,S 1,3=0 S1,4=1,S 1,5=−h 2,S 1,6=0(17.4.28 ) Similarly equation (17.4.24) yields S2,1=αkβk/2,S 2,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,S 2,6=S2,3 (17.4.29 ) while from equation (17.4.27) we find S3,1=0,S 3,2=0,S 3,3=−1 S3,4=0,S 3,5=0,S 3,6=1(17.4.30 ) Atx=0we have the boundary condition E3,1=/braceleftbiggy1,1,n−modd y2,1,n−meven(17.4.31 ) Recalltheconventionadoptedinthe solvderoutinethatforoneboundarycondition atk=1onlyS3,jcan be nonzero. Also, jtakes on the values 4 to 6 since the boundary condition involves only yk, not yk−1. Accordingly, the only nonzero values of S3,jatx=0are S3,4=1,n −modd S3,5=1,n −meven(17.4.32 ) 768 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Atx=1we 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 ),S 1,5=1,S 1,6=−y1,M 2(m+1 )(17.4.35) S2,4=1,S 2,5=0,S 2,6=0 (17.4.36) Herenowisthesampleprogramthatimplementstheabovealgorithm. Weneed a main program, sfroid, that calls the routine solvde, and we must supply the subroutine difeqcalled by solvde. For simplicity we choose an equally spaced 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=0does not involve y1ifn−mis even, we have to use the indexvfeature 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 y1andy2so that there is not a zero pivot element in s(i,j). The programprompts for values of mandn. It then computes an initial guess forybased on the Legendre function Pm n. It next prompts for c2, solves for y, promptsfor c2, solvesfor yusingthe previousvaluesas aninitial guess,andso on. PROGRAM sfroid INTEGER NE,M,NB,NCI,NCJ,NCK,NSI,NSJ,NYJ,NYKCOMMON /sfrcom/ x,h,mm,n,c2,anorm Communicates with difeq . PARAMETER (NE=3,M=41,NB=1,NCI=NE,NCJ=NE-NB+1,NCK=M+1,NSI=NE, * NSJ=2*NE+1,NYJ=NE,NYK=M) C USES plgndr,solvde Sample program using solvde . Computes eigenvalues of spheroidal harmonics Smn (x;c) form≥0andn≥m. In the program, mismm,c2isc2,a n d γof equation (17.4.20) isanorm . INTEGER i,itmax,k,mm,n,indexv(NE)REAL anorm,c2,conv,deriv,fac1,fac2,h,q1,slowc, * c(NCI,NCJ,NCK),s(NSI,NSJ),scalv(NE),x(M),y(NE,M),plgndr itmax=100conv=5.e-6 slowc=1. h=1./(M-1)c2=0.write(*,*)’ENTER M,N’ read(*,*)mm,n if(mod(n+mm,2).eq.1)then No interchanges necessary. indexv(1)=1 indexv(2)=2 indexv(3)=3 else Interchange y 1andy2. indexv(1)=2 indexv(2)=1 indexv(3)=3 endif anorm=1. Compute γ. 17.4AWorkedExample: SpheroidalHarmonics 769Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).if(mm.NE.0)then q1=n do11i=1,mm anorm=-.5*anorm*(n+i)*(q1/i)q1=q1-1. enddo 11 endif do12k=1,M-1 Initial guess. x(k)=(k-1)*h fac1=1.-x(k)**2 fac2=fac1**(-mm/2.)y(1,k)=plgndr(n,mm,x(k))*fac2 P m nfrom§6.8. deriv=-((n-mm+1)*plgndr(n+1,mm,x(k))-(n+1)* * x(k)*plgndr(n,mm,x(k)))/fac1 Derivative of Pm nfrom a recurrence re- lation. y(2,k)=mm*x(k)*y(1,k)/fac1+deriv*fac2 y(3,k)=n*(n+1)-mm*(mm+1) enddo 12 x(M)=1. Initial guess at x=1 done separately. y(1,M)=anormy(3,M)=n*(n+1)-mm*(mm+1) y(2,M)=(y(3,M)-c2)*y(1,M)/(2.*(mm+1.)) scalv(1)=abs(anorm)scalv(2)=max(abs(anorm),y(2,M)) scalv(3)=max(1.,y(3,M)) 1 continue write (*,*) ’ENTER C**2 OR 999 TO END’read (*,*) c2 if (c2.eq.999.) stop call solvde(itmax,conv,slowc,scalv,indexv,NE,NB,M,y,NYJ,NYK, * c,NCI,NCJ,NCK,s,NSI,NSJ) write (*,* )’M= ’,mm,’ N = ’,n, * ’ C**2 = ’,c2,’ LAMBDA = ’,y(3,1)+mm*(mm+1) goto 1 for another value of c 2. END SUBROUTINE difeq(k,k1,k2,jsf,is1,isf,indexv,ne,s,nsi,nsj,y,nyj,nyk) INTEGER is1,isf,jsf,k,k1,k2,ne,nsi,nsj,nyj,nyk,indexv(nyj),M REAL s(nsi,nsj),y(nyj,nyk)COMMON /sfrcom/ x,h,mm,n,c2,anormPARAMETER (M=41) Returns matrix s(i,j) forsolvde . INTEGER mm,nREAL anorm,c2,h,temp,temp2,x(M)if(k.eq.k1) then Boundary condition at first point. if(mod(n+mm,2).eq.1)then s(3,3+indexv(1))=1. Equation (17.4.32). s(3,3+indexv(2))=0. s(3,3+indexv(3))=0. s(3,jsf)=y(1,1) Equation (17.4.31). else s(3,3+indexv(1))=0. Equation (17.4.32). s(3,3+indexv(2))=1. s(3,3+indexv(3))=0.s(3,jsf)=y(2,1) Equation (17.4.31). endif else if(k.gt.k2) then Boundary conditions at last point. s(1,3+indexv(1))=-(y(3,M)-c2)/(2.*(mm+1.)) Equation (17.4.35). s(1,3+indexv(2))=1. s(1,3+indexv(3))=-y(1,M)/(2.*(mm+1.)) s(1,jsf)=y(2,M)-(y(3,M)-c2)*y(1,M)/(2.*(mm+1.)) Equation (17.4.33). s(2,3+indexv(1))=1. Equation (17.4.36). s(2,3+indexv(2))=0. 770 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).s(2,3+indexv(3))=0. s(2,jsf)=y(1,M)-anorm Equation (17.4.34). else Interior point. s(1,indexv(1))=-1. Equation (17.4.28). s(1,indexv(2))=-.5*h s(1,indexv(3))=0. s(1,3+indexv(1))=1.s(1,3+indexv(2))=-.5*hs(1,3+indexv(3))=0. temp=h/(1.-(x(k)+x(k-1))**2*.25) temp2=.5*(y(3,k)+y(3,k-1))-c2*.25*(x(k)+x(k-1))**2s(2,indexv(1))=temp*temp2*.5 Equation (17.4.29). s(2,indexv(2))=-1.-.5*temp*(mm+1.)*(x(k)+x(k-1)) s(2,indexv(3))=.25*temp*(y(1,k)+y(1,k-1)) s(2,3+indexv(1))=s(2,indexv(1))s(2,3+indexv(2))=2.+s(2,indexv(2)) s(2,3+indexv(3))=s(2,indexv(3)) s(3,indexv(1))=0. Equation (17.4.30). s(3,indexv(2))=0.s(3,indexv(3))=-1. s(3,3+indexv(1))=0. s(3,3+indexv(2))=0.s(3,3+indexv(3))=1. s(1,jsf)=y(1,k)-y(1,k-1)-.5*h*(y(2,k)+y(2,k-1)) Equation (17.4.23). s(2,jsf)=y(2,k)-y(2,k-1)-temp*((x(k)+x(k-1)) Equation (17.4.24). * *.5*(mm+1.)*(y(2,k)+y(2,k-1))-temp2** .5*(y(1,k)+y(1,k-1))) s(3,jsf)=y(3,k)-y(3,k-1) Equation (17.4.27). endifreturn END 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. SelectedOutputof sfroid mn c2λexact λsfroid 22 0 .16 .01427 6 .01427 1.06 .14095 6 .14095 4.06 .54250 6 .54253 25 1 .03 0 .4361 30 .4372 16.03 6 .9963 37 .0135 41 1 −1.0 131 .560 131 .554 Shooting Tosolvethesameproblemviashooting( §17.1),wesupplyasubroutine derivs thatimplementsequations(17.4.15)–(17.4.17). We will integratetheequationsover the range −1≤x≤0. We provide the subroutine loadwhich sets the eigenvalue y3to its current best estimate, v(1). It also sets the boundary values of y1and 17.4AWorkedExample: SpheroidalHarmonics 771Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).y2using 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 dxfrom the boundary to avoid having to evaluate y/prime 2right on the boundary. The subroutine scorefollows from equation (17.4.18). PROGRAM sphoot Sample program using shoot . Computes eigenvalues of spheroidal harmonics Smn (x;c)for m≥0andn≥m. Be sure that routine funcv fornewt is provided by shoot (§17.1). INTEGER i,m,n,nvar,N2 PARAMETER (N2=1) REAL c2,dx,gamma,q1,x1,x2,v(N2)LOGICAL check COMMON /sphcom/ c2,gamma,dx,m,n Communicates with load ,score ,a n d derivs . COMMON /caller/ x1,x2,nvar Communicates with shoot . C USES newt dx=1.e-4 Avoid evaluating derivatives exactly at x=−1. nvar=3 Number of equations. 1 write(*,*) ’input m,n,c-squared (999 to end)’ read(*,*) m,n,c2if (c2.eq.999.) stop if ((n.lt.m).or.(m.lt.0)) goto 1 gamma=1.0 Compute γof equation (17.4.20). q1=n do 11i=1,m gamma=-0.5*gamma*(n+i)*(q1/i)q1=q1-1.0 enddo 11 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 call newt(v,N2,check) Find vthat zeros function finscore . if(check)then write(*,*)’shoot failed; bad initial guess’ else write(*,’(1x,t6,a)’) ’mu(m,n)’ write(*,’(1x,f12.6)’) v(1)goto 1 endif END SUBROUTINE load(x1,v,y) INTEGER m,n REAL c2,dx,gamma,x1,y1,v(1),y(3)COMMON /sphcom/ c2,gamma,dx,m,n Supplies starting values for integration at x=−1+dx. y(3)=v(1) if(mod(n-m,2).eq.0)then y1=gamma else y1=-gamma endify(2)=-(y(3)-c2)*y1/(2*(m+1)) y(1)=y1+y(2)*dx returnEND SUBROUTINE score(x2,y,f) INTEGER m,nREAL c2,dx,gamma,x2,f(1),y(3) COMMON /sphcom/ c2,gamma,dx,m,n Tests whether boundary condition at x=0 is satisfied. if (mod(n-m,2).eq.0) then f(1)=y(2) 772 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).else f(1)=y(1) endif returnEND SUBROUTINE derivs(x,y,dydx) INTEGER m,nREAL c2,dx,gamma,x,dydx(3),y(3) COMMON /sphcom/ c2,gamma,dx,m,n 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 returnEND Shootingto a FittingPoint For variety we illustrate shootffrom§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 derivsis identical to the one for shoot. Now, however,there are two load routines. The routine load1forx=−1is essentially identical to loadabove. At x=1,load2sets the function value y1and the eigenvalue y3to their best current estimates, v2(1)andv2(2), respectively. If you quite sensibly make your initial guess of the eigenvaluethe same in the two intervals, then v1(1)will stay equal to v2(2)during the iteration. The subroutine scoresimply checks whether all three function values match at the fitting point. PROGRAM sphfpt Sample program using shootf . Computes eigenvalues of spheroidal harmonics Smn (x;c) form≥0andn≥m. Be sure that routine funcv fornewt is provided by shootf (§17.2). The routine derivs is the same as for sphoot . INTEGER i,m,n,nvar,nn2,N1,N2,NTOTREAL DXXPARAMETER (N1=2,N2=1,NTOT=N1+N2,DXX=1.e-4) REAL c2,dx,gamma,q1,x1,x2,xf,v1(N2),v2(N1),v(NTOT) LOGICAL checkCOMMON /sphcom/ c2,gamma,dx,m,n Communicates with load1 ,load2 ,score ,a n d derivs . COMMON /caller/ x1,x2,xf,nvar,nn2 Communicates with shootf . EQUIVALENCE (v1(1),v(1)),(v2(1),v(N2+1)) C USES newt nvar=NTOT Number of equations. nn2=N2dx=DXX Avoid evaluating derivatives exactly at x=±1. 1 write(*,*) ’input m,n,c-squared (999 to end)’ read(*,*) m,n,c2 if (c2.eq.999.) stopif ((n.lt.m).or.(m.lt.0)) goto 1 gamma=1.0 Compute γof equation (17.4.20). q1=ndo 11i=1,m gamma=-0.5*gamma*(n+i)*(q1/i) q1=q1-1.0 enddo 11 v1(1)=n*(n+1)-m*(m+1)+c2/2.0 Initial guess for eigenvalue and function value. v2(2)=v1(1) 17.4A WorkedExample: SpheroidalHarmonics 773Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).v2(1)=gamma*(1.-(v2(2)-c2)*dx/(2*(m+1))) x1=-1.0+dx Set range of integration. x2=1.0-dx xf=0. Fitting point. call newt(v,NTOT,check) Find vthat zeros function finscore . if(check)then write(*,*)’shootf failed; bad initial guess’ else write(*,’(1x,t6,a)’) ’mu(m,n)’ write(*,’(1x,f12.6)’) v1(1) goto 1 endifEND SUBROUTINE load1(x1,v1,y) INTEGER m,n REAL c2,dx,gamma,x1,y1,v1(1),y(3) COMMON /sphcom/ c2,gamma,dx,m,n Supplies starting values for integration at x=−1+dx. y(3)=v1(1) if(mod(n-m,2).eq.0)then y1=gamma else y1=-gamma endify(2)=-(y(3)-c2)*y1/(2*(m+1))y(1)=y1+y(2)*dx return END SUBROUTINE load2(x2,v2,y) INTEGER m,nREAL c2,dx,gamma,x2,v2(2),y(3)COMMON /sphcom/ c2,gamma,dx,m,n Supplies starting values for integration at x=1−dx. y(3)=v2(2)y(1)=v2(1)y(2)=(y(3)-c2)*y(1)/(2*(m+1)) return END SUBROUTINE score(xf,y,f) INTEGER i,m,nREAL c2,gamma,dx,xf,f(3),y(3)COMMON /sphcom/ c2,gamma,dx,m,n Tests whether solutions match at fitting point x=0. do 12i=1,3 f(i)=y(i) enddo 12 return END CITED REFERENCES AND FURTHER READING: Flammer, C. 1957, Spheroidal Wave Functions (Stanford, CA: Stanford University Press). [1] Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §21. [2] Morse,P.M.,andFeshbach,H.1953, MethodsofTheoreticalPhysics ,PartII(NewYork:McGraw- Hill), pp. 1502ff. [3] 774 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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 thenumber of points usually cures the problem. If we know ahead of time where our solutions will be rapidly varying, we can put more gridpointsthereandlesselsewhere. Alternatively,wecansolvetheproblemfirstonauniformmesh and then examine the solution to see where we should add more points. We then repeatthe 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 improvesthe accuracy of the relaxation method, but also (as we will see in the next section) allowsinternal singularities to be handled in quite a neat way. Here we learn how to accomplishthe automatic allocation. We want to focus attention on the independent variable x, and consider two alternative reparametrizations of it. The first, we term q; this is just the coordinate corresponding to the mesh points themselves, so that q=1atk=1,q=2atk=2, and so on. Between any two meshpointswehave ∆q=1. Inthechangeofindependent variableintheODEsfrom xtoq, dy dx=g (17.5.1 ) becomes dy dq=gdx dq(17.5.2 ) In terms of q, equation (17.5.2) as an FDE might be written yk−yk−1−1 2/bracketleftBigg/parenleftBigg gdx dq/parenrightBigg k+/parenleftBigg gdx dq/parenrightBigg k−1/bracketrightBigg =0 ( 17.5.3 ) or some related version. Note that dx/dqshould accompany g. The transformation between xandqdepends only on the Jacobian dx/dq. Its reciprocal dq/dxis 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/dxto be larger where yis changing rapidly, or near to the boundaries, or both. In fact, we can probably make up a formula for what we would like dq/dxtobeproportionalto. Theproblemisthatwedonotknowtheproportionalityconstant. That is, the formula that we might invent would not have the correct integral over the wholerangeof xsoastomake qvaryfrom 1toM,accordingtoitsdefinition. Tosolvethisproblem we introduce a second reparametrization Q(q), where Qis a new independent variable. The relation between Qandqis taken to be linear, so that a mesh spacing formula for dQ/dx differs only in its unknown proportionality constant. A linear relation implies d 2Q dq2=0 ( 17.5.4 ) or, expressed in the usual manner as coupled first-order equations, dQ(x) dq=ψdψ dq=0 ( 17.5.5 ) 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 φ(x)=dQ dx=dQ dqdq dx(17.5.6 )