c17-2 (779602)
Текст из файла
760Chapter 17.Two Point Boundary Value Problems17.2 Shooting to a Fitting Pointyi (x1 ) = yi (x1 ; V(1)1 , . . . , V(1)n2 )i = 1, . . . , N(17.2.1)Likewise we can define an n1 -vector V(2) of starting parameters at x2 , and aprescription load2(x2,v2,y) for mapping V(2) into a y that satisfies the boundaryconditions at x2 ,yi (x2 ) = yi (x2 ; V(2)1 , .
. . , V(2)n1 )i = 1, . . . , N(17.2.2)We thus have a total of N freely adjustable parameters in the combination ofV(1) and V(2) . The N conditions that must be satisfied are that there be agreementin N components of y at xf between the values obtained integrating from one sideand from the other,yi (xf ; V(1) ) = yi (xf ; V(2) )i = 1, . . .
, N(17.2.3)In some problems, the N matching conditions can be better described (physically,mathematically, or numerically) by using N different functions Fi , i = 1 . . . N , eachpossibly depending on the N components yi . In those cases, (17.2.3) is replaced byFi [y(xf ; V(1))] = Fi [y(xf ; V(2))]i = 1, . . . , N(17.2.4)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).The shooting method described in §17.1 tacitly assumed that the “shots” wouldbe able to traverse the entire domain of integration, even at the early stages ofconvergence to a correct solution. In some problems it can happen that, for verywrong starting conditions, an initial solution can’t even get from x1 to x2 withoutencountering some incalculable, or catastrophic, result.
For example, the argumentof a square root might go negative, causing the numerical code to crash. Simpleshooting would be stymied.A different, but related, case is where the endpoints are both singular pointsof the set of ODEs. One frequently needs to use special methods to integrate nearthe singular points, analytic asymptotic expansions, for example.
In such cases it isfeasible to integrate in the direction away from a singular point, using the specialmethod to get through the first little bit and then reading off “initial” values forfurther numerical integration. However it is usually not feasible to integrate intoa singular point, if only because one has not usually expended the same analyticeffort to obtain expansions of “wrong” solutions near the singular point (those notsatisfying the desired boundary condition).The solution to the above mentioned difficulties is shooting to a fitting point.Instead of integrating from x1 to x2 , we integrate first from x1 to some point xf thatis between x1 and x2 ; and second from x2 (in the opposite direction) to xf .If (as before) the number of boundary conditions imposed at x1 is n1 , and thenumber imposed at x2 is n2 , then there are n2 freely specifiable starting values atx1 and n1 freely specifiable starting values at x2 . (If you are confused by this, goback to §17.1.) We can therefore define an n2 -vector V(1) of starting parametersat x1 , and a prescription load1(x1,v1,y) for mapping V(1) into a y that satisfiesthe boundary conditions at x1 ,17.2 Shooting to a Fitting Point761#include "nrutil.h"#define EPS 1.0e-6extern int nn2,nvar;extern float x1,x2,xf;Variables that you must define and set in your main program.int kmax,kount;float *xp,**yp,dxsav;Communicates with odeint.void shootf(int n, float v[], float f[])Routine for use with newt to solve a two point boundary value problem for nvar coupledODEs by shooting from x1 and x2 to a fitting point xf.
Initial values for the nvar ODEs atx1 (x2) are generated from the n2 (n1) coefficients v1 (v2), using the user-supplied routineload1 (load2). The coefficients v1 and v2 should be stored in a single array v[1..n1+n2]in the main program by statements of the form v1=v; and v2 = &v[n2];. The input parameter n = n1 + n2 = nvar. The routine integrates the ODEs to xf using the Runge-Kuttamethod with tolerance EPS, initial stepsize h1, and minimum stepsize hmin. At xf it calls theuser-supplied routine score to evaluate the nvar functions f1 and f2 that ought to matchat xf. The differences f are returned on output.
newt uses a globally convergent Newton’smethod to adjust the values of v until the functions f are zero. The user-supplied routinederivs(x,y,dydx) supplies derivative information to the ODE integrator (see Chapter 16).The first set of global variables above receives its values from the main program so that shootcan have the syntax required for it to be the argument vecfunc of newt. Set nn2 = n2 inthe main program.{void derivs(float x, float y[], float dydx[]);void load1(float x1, float v1[], float y[]);void load2(float x2, float v2[], 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,float [], 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 i,nbad,nok;float h1,hmin=0.0,*f1,*f2,*y;f1=vector(1,nvar);f2=vector(1,nvar);y=vector(1,nvar);kmax=0;h1=(x2-x1)/100.0;load1(x1,v,y);Path from x1 to xf with best trial values v1.odeint(y,nvar,x1,xf,EPS,h1,hmin,&nok,&nbad,derivs,rkqs);score(xf,y,f1);load2(x2,&v[nn2],y);Path from x2 to xf with best trial values v2.odeint(y,nvar,x2,xf,EPS,h1,hmin,&nok,&nbad,derivs,rkqs);score(xf,y,f2);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 the program below, the user-supplied function score(xf,y,f) is supposedto map an input N -vector y into an output N -vector F.
In most cases, you candummy this function as the identity mapping.Shooting to a fitting point uses globally convergent Newton-Raphson exactlyas in §17.1. Comparing closely with the routine shoot of the previous section, youshould have no difficulty in understanding the following routine shootf.
The maindifferences in use are that you have to supply both load1 and load2. Also, in thecalling program you must supply initial guesses for v1[1..n2] and v2[1..n1].Once again a sample program illustrating shooting to a fitting point is given in §17.4.762Chapter 17.Two Point Boundary Value Problemsfor (i=1;i<=n;i++) f[i]=f1[i]-f2[i];free_vector(y,1,nvar);free_vector(f2,1,nvar);free_vector(f1,1,nvar);}CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America).Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA:Blaisdell).Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§§7.3.5–7.3.6.
[1]17.3 Relaxation MethodsIn relaxation methods we replace ODEs by approximate finite-difference equations(FDEs) on a grid or mesh of points that spans the domain of interest. As a typical example,we could replace a general first-order differential equationdy= g(x, y)dxwith an algebraic equation relating function values at two points k, k − 1:yk − yk−1 − (xk − xk−1 ) g 12 (xk + xk−1 ), 12 (yk + yk−1 ) = 0(17.3.1)(17.3.2)The form of the FDE in (17.3.2) illustrates the idea, but not uniquely: There are manyways to turn the ODE into an FDE. When the problem involves N coupled first-order ODEsrepresented by FDEs on a mesh of M points, a solution consists of values for N dependentfunctions given at each of the M mesh points, or N × M variables in all. The relaxationmethod determines the solution by starting with a guess and improving it, iteratively.
As theiterations improve the solution, the result is said to relax to the true solution.While several iteration schemes are possible, for most problems our old standby, multidimensional Newton’s method, works well. The method produces a matrix equation thatmust be solved, but the matrix takes a special, “block diagonal” form, that allows it to beinverted far more economically both in time and storage than would be possible for a generalmatrix of size (M N ) × (M N ).
Since M N can easily be several thousand, this is crucialfor the feasibility of the method.Our implementation couples at most pairs of points, as in equation(17.3.2). More points can be coupled, but then the method becomes more complex.We will provide enough background so that you can write a more general scheme if youhave the patience to do so.Let us develop a general set of algebraic equations that represent the ODEs by FDEs. TheODE problem is exactly identical to that expressed in equations (17.0.1)–(17.0.3) where wehad N coupled first-order equations that satisfy n1 boundary conditions at x1 and n2 = N −n1boundary conditions at x2 .
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
