c17-4 (779604), страница 2
Текст из файла (страница 2)
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit 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).Substituting this expansion in equation (17.4.4) and letting x → 1, we find that17.4 A Worked Example: Spheroidal Harmonics775The boundary condition at x = 0 in this notation isy1 = 0,n − m oddy2 = 0,n − m even(17.4.18)y2 =y3 − c2y12(m + 1)(17.4.19)y1 = lim (1 − x2 )−m/2 Pnm (x) =x→1(−1)m (n + m)!≡γ2m m!(n − m)!(17.4.20)We are now ready to illustrate the use of the methods of previous sectionson this problem.RelaxationIf we just want a few isolated values of λ or S, shooting is probably the quickestmethod.
However, if we want values for a large sequence of values of c, relaxationis better. Relaxation rewards a good initial guess with rapid convergence, and theprevious 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 atotal of M mesh points, we haveh=1M −1xk = (k − 1)h,(17.4.21)k = 1, 2, . . .
, M(17.4.22)At interior points k = 2, 3, . . ., M , equation (17.4.15) givesE1,k = y1,k − y1,k−1 −h(y2,k + y2,k−1 )2(17.4.23)Equation (17.4.16) givesE2,k = y2,k − y2,k−1 − βk(y1,k + y1,k−1 )(xk + xk−1 )(m + 1)(y2,k + y2,k−1 )− αk×22(17.4.24)wherey3,k + y3,k−1 c2 (xk + xk−1)2−24hβk =1 − 14 (xk + xk−1 )2αk =(17.4.25)(17.4.26)Finally, equation (17.4.17) givesE3,k = y3,k − y3,k−1(17.4.27)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit 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).At x = 1 we have two conditions:776Chapter 17.Two Point Boundary Value ProblemsNow recall that the matrix of partial derivatives Si,j of equation (17.3.8) isdefined so that i labels the equation and j the variable. In our case, j runs from 1 to3 for yj at k − 1 and from 4 to 6 for yj at k.
Thus equation (17.4.23) givesS1,1 = −1,S1,3 = 0(17.4.28)S1,6 = 0Similarly equation (17.4.24) yieldsS2,1 = αk βk /2,S2,2 = −1 − βk (xk + xk−1)(m + 1)/2,S2,3 = βk (y1,k + y1,k−1)/4S2,5 = 2 + S2,2 ,S2,4 = S2,1 ,S2,6 = S2,3(17.4.29)while from equation (17.4.27) we findS3,1 = 0,S3,4 = 0,S3,2 = 0,S3,5 = 0,S3,3 = −1S3,6 = 1(17.4.30)At x = 0 we have the boundary conditiony1,1 , n − m oddE3,1 =y2,1 , n − m even(17.4.31)Recall the convention adopted in the solvde routine that for one boundary conditionat k = 1 only S3,j can be nonzero. Also, j takes on the values 4 to 6 since theboundary condition involves only yk , not yk−1 . Accordingly, the only nonzerovalues of S3,j at x = 0 areS3,4 = 1,n − m oddS3,5 = 1,n − m even(17.4.32)At x = 1 we havey3,M − c2y1,M2(m + 1)−γE1,M +1 = y2,M −(17.4.33)E2,M +1 = y1,M(17.4.34)ThusS1,4 = −S2,4 = 1,y3,M − c2,2(m + 1)S1,5 = 1,S1,6 = −S2,5 = 0,S2,6 = 0y1,M2(m + 1)(17.4.35)(17.4.36)Here now is the sample program that implements the above algorithm.
Weneed a main program, sfroid, that calls the routine solvde, and we must supplythe function difeq called by solvde. For simplicity we choose an equally spacedSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit 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).S1,4 = 1,hS1,2 = − ,2hS1,5 = − ,217.4 A Worked Example: Spheroidal Harmonics777#include <stdio.h>#include <math.h>#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) form ≥ 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) {indexv[1]=1;indexv[2]=2;indexv[3]=3;} else {indexv[1]=2;indexv[2]=1;indexv[3]=3;}anorm=1.0;if (mm) {q1=n;No interchanges necessary.Interchange y1 and y2 .Compute γ.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit 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).mesh of m = 41 points, that is, h = .025. As we shall see, this gives good accuracyfor the eigenvalues up to moderate values of n − m.Since the boundary condition at x = 0 does not involve y1 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 − mis even, we need to interchange the columns for y1 and y2 so that there is not azero pivot element in s[i][j].The program prompts for values of m and n. It then computes an initial guessfor y based on the Legendre function Pnm . 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.778Chapter 17.Two Point Boundary Value Problems}extern int mm,n,mpt;extern float h,c2,anorm,x[];Defined in sfroid.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).Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit 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).for (i=1;i<=mm;i++) anorm = -0.5*anorm*(n+i)*(q1--/i);}for (k=1;k<=(M-1);k++) {Initial guess.x[k]=(k-1)*h;fac1=1.0-x[k]*x[k];fac2=exp((-mm/2.0)*log(fac1));y[1][k]=plgndr(n,mm,x[k])*fac2;Pnm from §6.8.deriv = -((n-mm+1)*plgndr(n+1,mm,x[k])Derivative of Pnm from a recur(n+1)*x[k]*plgndr(n,mm,x[k]))/fac1;rence relation.y[2][k]=mm*x[k]*y[1][k]/fac1+deriv*fac2;y[3][k]=n*(n+1)-mm*(mm+1);}x[M]=1.0;Initial guess at x = 1 done sepy[1][M]=anorm;arately.y[3][M]=n*(n+1)-mm*(mm+1);y[2][M]=(y[3][M]-c2)*y[1][M]/(2.0*(mm+1.0));scalv[1]=fabs(anorm);scalv[2]=(y[2][M] > 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 .17.4 A Worked Example: Spheroidal Harmonics(17.4.23).(17.4.24).}You can run the program and check it against values of λmn (c) given inthe tables at the back of Flammer’s book [1] or in Table 21.1 of Abramowitz andStegun [2].
Typically it converges in about 3 iterations. The table below givesa few comparisons.mSelected Output of sfroidnc2λexactλsfroid22254110.11.04.01.016.0−1.06.014276.140956.5425030.436136.9963131.5606.014276.140956.5425330.437237.0135131.554ShootingTo solve the same problem via shooting (§17.1), we supply a function derivsthat implements equations (17.4.15)–(17.4.17). We will integrate the equations overthe range −1 ≤ x ≤ 0. We provide the function load which sets the eigenvaluey3 to its current best estimate, v[1]. It also sets the boundary values of y1 andy2 using equations (17.4.20) and (17.4.19) (with a minus sign corresponding toSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit 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).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]);s[2][jsf]=y[2][k]-y[2][k-1]-temp*((x[k]+x[k-1])*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).}779780Chapter 17.Two Point Boundary Value Problemsx = −1).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
