c17-1 (779601)
Текст из файла
75717.1 The Shooting Method17.1 The Shooting Methodyi (x1 ) = yi (x1 ; V1 , . . . , Vn2 )i = 1, . . . , N(17.1.1)Below, the function that implements (17.1.1) will be called load.Notice that the components of V might be exactly the values of certain “free”components of y, with the other components of y determined by the boundaryconditions. Alternatively, the components of V might parametrize the solutions thatsatisfy the starting boundary conditions in some other convenient way. Boundaryconditions often impose algebraic relations among the yi , rather than specific valuesfor each of them. Using some auxiliary set of parameters often makes it easier to“solve” the boundary relations for a consistent set of yi ’s.
It makes no differencewhich way you go, as long as your vector space of V’s generates (through 17.1.1)all allowed starting vectors y.Given a particular V, a particular y(x1 ) is thus generated. It can then be turnedinto a y(x2 ) by integrating the ODEs to x2 as an initial value problem (e.g., usingChapter 16’s odeint). Now, at x2 , let us define a discrepancy vector F, also ofdimension n2 , whose components measure how far we are from satisfying the n2boundary conditions at x2 (17.0.3). Simplest of all is just to use the right-handsides of (17.0.3),Fk = B2k (x2 , y)k = 1, . . . , n2(17.1.2)As in the case of V, however, you can use any other convenient parametrization,as long as your space of F’s spans the space of possible discrepancies from thedesired boundary conditions, with all components of F equal to zero if and only ifthe boundary conditions at x2 are satisfied.
Below, you will be asked to supply auser-written function score which uses (17.0.3) to convert an N -vector of endingvalues y(x2 ) into an n2 -vector of discrepancies F.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 this section we discuss “pure” shooting, where the integration proceeds fromx1 to x2 , and we try to match boundary conditions at the end of the integration. Inthe next section, we describe shooting to an intermediate fitting point, where thesolution to the equations and boundary conditions is found by launching “shots”from both sides of the interval and trying to match continuity conditions at someintermediate point.Our implementation of the shooting method exactly implements multidimensional, globally convergent Newton-Raphson (§9.7). It seeks to zero n2 functionsof n2 variables.
The functions are obtained by integrating N differential equationsfrom x1 to x2 . Let us see how this works:At the starting point x1 there are N starting values yi to be specified, butsubject to n1 conditions. Therefore there are n2 = N − n1 freely specifiable startingvalues. Let us imagine that these freely specifiable values are the components of avector V that lives in a vector space of dimension n2 .
Then you, the user, knowingthe functional form of the boundary conditions (17.0.2), can write a function thatgenerates a complete set of N starting values y, satisfying the boundary conditionsat x1 , from an arbitrary vector value of V in which there are no restrictions on the n2component values. In other words, (17.0.2) converts to a prescription758Chapter 17.Two Point Boundary Value ProblemsNow, as far as Newton-Raphson is concerned, we are nearly in business.
Wewant to find a vector value of V that zeros the vector value of F. We do thisby invoking the globally convergent Newton’s method implemented in the routinenewt of §9.7. Recall that the heart of Newton’s method involves solving the setof n2 linear equations(17.1.3)and then adding the correction back,Vnew = Vold + δV(17.1.4)In (17.1.3), the Jacobian matrix J has components given byJij =∂Fi∂Vj(17.1.5)It is not feasible to compute these partial derivatives analytically.
Rather, eachrequires a separate integration of the N ODEs, followed by the evaluation ofFi (V1 , . . . , Vj + ∆Vj , . . .) − Fi (V1 , . . . , Vj , . . .)∂Fi≈∂Vj∆Vj(17.1.6)This is done automatically for you in the routine fdjac that comes with newt. Theonly input to newt that you have to provide is the routine vecfunc that calculatesF by integrating the ODEs. Here is the appropriate routine, called shoot, that isto be passed as the actual argument in newt:#include "nrutil.h"#define EPS 1.0e-6extern int nvar;extern float x1,x2;Variables that you must define and set in your main program.int kmax,kount;float *xp,**yp,dxsav;Communicates with odeint.void shoot(int n, float v[], float f[])Routine for use with newt to solve a two point boundary value problem for nvar coupled ODEsby shooting from x1 to x2.
Initial values for the nvar ODEs at x1 are generated from the n2input coefficients v[1..n2], using the user-supplied routine load. The routine integrates theODEs to x2 using the Runge-Kutta method with tolerance EPS, initial stepsize h1, and minimumstepsize hmin. At x2 it calls the user-supplied routine score to evaluate the n2 functionsf[1..n2] that ought to be zero to satisfy the boundary conditions at x2. The functions fare returned on output. newt uses a globally convergent Newton’s method to adjust the valuesof v until the functions f are zero.
The user-supplied routine derivs(x,y,dydx) suppliesderivative information to the ODE integrator (see Chapter 16). The first set of global variablesabove receives its values from the main program so that shoot can have the syntax requiredfor it to be the argument vecfunc of newt.{void derivs(float x, float y[], float dydx[]);void load(float x1, float v[], float y[]);void odeint(float ystart[], int nvar, float x1, float x2,float eps, float h1, float hmin, int *nok, int *nbad,void (*derivs)(float, float [], float []),void (*rkqs)(float [], float [], int, float *, float, float,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).J · δV = −F17.1 The Shooting Method759float [], float *, float *, void (*)(float, float [], float [])));void rkqs(float y[], float dydx[], int n, float *x,float htry, float eps, float yscal[], float *hdid, float *hnext,void (*derivs)(float, float [], float []));void score(float xf, float y[], float f[]);int nbad,nok;float h1,hmin=0.0,*y;}For some problems the initial stepsize ∆V might depend sensitively upon theinitial conditions.
It is straightforward to alter load to include a suggested stepsizeh1 as another output variable and feed it to fdjac via a global variable.A complete cycle of the shooting method thus requires n2 + 1 integrations ofthe N coupled ODEs: one integration to evaluate the current degree of mismatch,and n2 for the partial derivatives. Each new cycle requires a new round of n2 + 1integrations. This illustrates the enormous extra effort involved in solving two pointboundary value problems compared with intial value problems.If the differential equations are linear, then only one complete cycle is required,since (17.1.3)–(17.1.4) should take us right to the solution.
A second round can beuseful, however, in mopping up some (never all) of the roundoff error.As given here, shoot uses the quality controlled Runge-Kutta method of §16.2to integrate the ODEs, but any of the other methods of Chapter 16 could just aswell be used.You, the user, must supply shoot with: (i) a function load(x1,v,y) whichcalculates the n-vector y[1..n] (satisfying the starting boundary conditions, ofcourse), given the freely specifiable variables of v[1..n2] at the initial point x1;(ii) a function score(x2,y,f) which calculates the discrepancy vector f[1..n2]of the ending boundary conditions, given the vector y[1..n] at the endpoint x2;(iii) a starting vector v[1..n2]; (iv) a function derivs for the ODE integration; andother obvious parameters as described in the header comment above.In §17.4 we give a sample program illustrating how to use shoot.CITED REFERENCES AND FURTHER READING:Acton, F.S.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.















