c17-4 (779604), страница 3
Текст из файла (страница 3)
Note that the boundary condition is actually applied a distance dx fromthe boundary to avoid having to evaluate y20 right on the boundary. The functionscore follows from equation (17.4.18).int m,n;float c2,dx,gmma;Communicates with load, score, and derivs.int nvar;float x1,x2;Communicates with shoot.int main(void) /* Program sphoot */Sample program using shoot. Computes eigenvalues of spheroidal harmonics Smn (x; c) form ≥ 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]);}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).#include <stdio.h>#include "nrutil.h"#define N2 117.4 A Worked Example: Spheroidal Harmonics781Shooting to a Fitting PointFor variety we illustrate shootf from §17.2 by integrating over the wholerange −1 + dx ≤ x ≤ 1 − dx, with the fitting point chosen to be at x = 0.
Theroutine derivs is identical to the one for shoot. Now, however, there are two loadroutines. The routine load1 for x = −1 is essentially identical to load above. Atx = 1, load2 sets the function value y1 and the eigenvalue y3 to their best currentestimates, v2[1] and v2[2], respectively. If you quite sensibly make your initialguess of the eigenvalue the same in the two intervals, then v1[1] will stay equalto v2[2] during the iteration. The function score simply checks whether all threefunction values match at the fitting point.#include <stdio.h>#include <math.h>#include "nrutil.h"#define N1 2#define N2 1#define NTOT (N1+N2)#define DXX 1.0e-4int m,n;float c2,dx,gmma;Communicates with load1, load2, score,and derivs.int nn2,nvar;float x1,x2,xf;Communicates with shootf.int main(void) /* Program sphfpt */Sample program using shootf.
Computes eigenvalues of spheroidal harmonics Smn (x; c) form ≥ 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;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).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;}782Chapter 17.Two Point Boundary Value Problems}free_vector(v,1,NTOT);return 0;}void load1(float x1, float v1[], float y[])Supplies starting values for integration at x = −1 + dx.{float y1 = (n-m & 1 ? -gmma : gmma);y[3]=v1[1];y[2] = -(y[3]-c2)*y1/(2*(m+1));y[1]=y1+y[2]*dx;}void load2(float x2, float v2[], float y[])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));}void score(float xf, float y[], float f[])Tests whether solutions match at fitting point x = 0.{int i;for (i=1;i<=3;i++) f[i]=y[i];}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 Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §21. [2]Morse, P.M., and Feshbach, H. 1953, Methods of Theoretical Physics, Part II (New York: McGrawHill), pp. 1502ff. [3]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<=m;i++) gmma *= -0.5*(n+i)*(q1--/i);v1[1]=n*(n+1)-m*(m+1)+c2/2.0;Initial guess for eigenvalue and function value.v2[2]=v1[1];v2[1]=gmma*(1.0-(v2[2]-c2)*dx/(2*(m+1)));x1 = -1.0+dx;Set range of integration.x2=1.0-dx;xf=0.0;Fitting point.newt(v,NTOT,&check,shootf);Find v that zeros function f in score.if (check) {printf("shootf failed; bad initial guess\n");} else {printf("\tmu(m,n)\n");printf("%12.6f\n",v[1]);}78317.5 Automated Allocation of Mesh Points17.5 Automated Allocation of Mesh Pointsdy=gdx(17.5.1)becomesdydx=gdqdqIn terms of q, equation (17.5.2) as an FDE might be written"!!dxdx1yk − yk−1 − 2g+ gdqdqk(17.5.2)#=0(17.5.3)k−1or some related version.
Note that dx/dq should accompany g. The transformation betweenx and q depends only on the Jacobian dx/dq. Its reciprocal dq/dx is proportional to thedensity 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 theboundaries, or both. In fact, we can probably make up a formula for what we would likedq/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 wholerange of x so as to make q vary from 1 to M , according to its definition. To solve this problemwe introduce a second reparametrization Q(q), where Q is a new independent variable. Therelation between Q and q is taken to be linear, so that a mesh spacing formula for dQ/dxdiffers only in its unknown proportionality constant. A linear relation impliesd2 Q=0dq2or, expressed in the usual manner as coupled first-order equations,(17.5.4)dQ(x)dψ=ψ=0(17.5.5)dqdqwhere ψ is a new intermediate variable. We add these two equations to the set of ODEsbeing solved.Completing the prescription, we add a third ODE that is just our desired mesh-densityfunction, namelyφ(x) =dQdQ dq=dxdq dx(17.5.6)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).In relaxation problems, you have to choose values for the independent variable at themesh points.
This is called allocating the grid or mesh. The usual procedure is to picka 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 moregrid points there and less elsewhere. Alternatively, we can solve the problem first on a uniformmesh 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 sucha 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 alternativereparametrizations of it. The first, we term q; this is just the coordinate corresponding to themesh points themselves, so that q = 1 at k = 1, q = 2 at k = 2, and so on.
Between any twomesh points we have ∆q = 1. In the change of independent variable in the ODEs from x to q,.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
