Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 83
Текст из файла (страница 83)
You should weigh timeagainst memory and convenience carefully.Of course neither of the above routines should be used for the trivial cases offinding the largest, or smallest, element in an array. Those cases, you code by handas simple for loops. There are also good ways to code the case where k is modest incomparison to N , so that extra memory of order k is not burdensome. An exampleis to use the method of Heapsort (§8.3) to make a single pass through an array oflength N while saving the m largest elements. The advantage of the heap structureis that only log m, rather than m, comparisons are required every time a new√elementis added to the candidate list. This becomes a real savings when m > O( N), butit never hurts otherwise and is easy to code. The following program gives the idea.void hpsel(unsigned long m, unsigned long n, float arr[], float heap[])Returns in heap[1..m] the largest m elements of the array arr[1..n], with heap[1] guaranteed to be the the mth largest element.
The array arr is not altered. For efficiency, this routineshould be used only when m n.{void sort(unsigned long n, float arr[]);void nrerror(char error_text[]);unsigned long i,j,k;float swap;if (m > n/2 || m < 1) nrerror("probable misuse of hpsel");for (i=1;i<=m;i++) heap[i]=arr[i];sort(m,heap);Create initial heap by overkill! We assume m n.for (i=m+1;i<=n;i++) {For each remaining element...if (arr[i] > heap[1]) {Put it on the heap?heap[1]=arr[i];for (j=1;;) {Sift down.8.6 Determination of Equivalence Classes345k=j << 1;if (k > m) break;if (k != m && heap[k] > heap[k+1]) k++;if (heap[j] <= heap[k]) break;swap=heap[k];heap[k]=heap[j];heap[j]=swap;j=k;}}}}CITED REFERENCES AND FURTHER READING:Sedgewick, R. 1988, Algorithms, 2nd ed. (Reading, MA: Addison-Wesley), pp. 126ff.
[1]Knuth, D.E. 1973, Sorting and Searching, vol. 3 of The Art of Computer Programming (Reading,MA: Addison-Wesley).8.6 Determination of Equivalence ClassesA number of techniques for sorting and searching relate to data structures whose detailsare beyond the scope of this book, for example, trees, linked lists, etc. These structures andtheir manipulations are the bread and butter of computer science, as distinct from numericalanalysis, and there is no shortage of books on the subject.In working with experimental data, we have found that one particular such manipulation,namely the determination of equivalence classes, arises sufficiently often to justify inclusionhere.The problem is this: There are N “elements” (or “data points” or whatever), numbered1, . . .
, N . You are given pairwise information about whether elements are in the sameequivalence class of “sameness,” by whatever criterion happens to be of interest. Forexample, you may have a list of facts like: “Element 3 and element 7 are in the same class;element 19 and element 4 are in the same class; element 7 and element 12 are in the sameclass, . . . .” Alternatively, you may have a procedure, given the numbers of two elementsj and k, for deciding whether they are in the same class or different classes. (Recall thatan equivalence relation can be anything satisfying the RST properties: reflexive, symmetric,transitive. This is compatible with any intuitive definition of “sameness.”)The desired output is an assignment to each of the N elements of an equivalence classnumber, such that two elements are in the same class if and only if they are assigned thesame class number.Efficient algorithms work like this: Let F (j) be the class or “family” number of elementj.
Start off with each element in its own family, so that F (j) = j. The array F (j) can beinterpreted as a tree structure, where F (j) denotes the parent of j. If we arrange for each familyto be its own tree, disjoint from all the other “family trees,” then we can label each family(equivalence class) by its most senior great-great-. . .grandparent. The detailed topology ofthe tree doesn’t matter at all, as long as we graft each related element onto it somewhere.Therefore, we process each elemental datum “j is equivalent to k” by (i) tracking jup to its highest ancestor, (ii) tracking k up to its highest ancestor, (iii) giving j to k as anew parent, or vice versa (it makes no difference). After processing all the relations, we gothrough all the elements j and reset their F (j)’s to their highest possible ancestors, whichthen label the equivalence classes.The following routine, based on Knuth [1], assumes that there are m elemental piecesof information, stored in two arrays of length m, lista,listb, the interpretation beingthat lista[j] and listb[j], j=1...m, are the numbers of two elements which (we arethus told) are related.346Chapter 8.Sortingvoid eclass(int nf[], int n, int lista[], int listb[], int m)Given m equivalences between pairs of n individual elements in the form of the input arrayslista[1..m] and listb[1..m], this routine returns in nf[1..n] the number of the equivalence class of each of the n elements, integers between 1 and n (not all such integers used).{int l,k,j;for (k=1;k<=n;k++) nf[k]=k;Initialize each element its own class.for (l=1;l<=m;l++) {For each piece of input information...j=lista[l];while (nf[j] != j) j=nf[j];Track first element up to its ancestor.k=listb[l];while (nf[k] != k) k=nf[k];Track second element up to its ancestor.if (j != k) nf[j]=k;If they are not already related, make them}so.for (j=1;j<=n;j++)Final sweep up to highest ancestors.while (nf[j] != nf[nf[j]]) nf[j]=nf[nf[j]];}Alternatively, we may be able to construct a function equiv(j,k) that returns a nonzero(true) value if elements j and k are related, or a zero (false) value if they are not.
Then wewant to loop over all pairs of elements to get the complete picture. D. Eardley has deviseda clever way of doing this while simultaneously sweeping the tree up to high ancestors in amanner that keeps it current and obviates most of the final sweep phase:void eclazz(int nf[], int n, int (*equiv)(int, int))Given a user-supplied boolean function equiv which tells whether a pair of elements, each inthe range 1...n, are related, return in nf[1..n] equivalence class numbers for each element.{int kk,jj;}nf[1]=1;for (jj=2;jj<=n;jj++) {Loop over first element of all pairs.nf[jj]=jj;for (kk=1;kk<=(jj-1);kk++) {Loop over second element of all pairs.nf[kk]=nf[nf[kk]];Sweep it up this much.if ((*equiv)(jj,kk)) nf[nf[nf[kk]]]=jj;Good exercise for the reader to figure out why this much ancestry is necessary!}}for (jj=1;jj<=n;jj++) nf[jj]=nf[nf[jj]];Only this much sweeping is neededfinally.CITED REFERENCES AND FURTHER READING:Knuth, D.E.
1968, Fundamental Algorithms, vol. 1 of The Art of Computer Programming (Reading,MA: Addison-Wesley), §2.3.3. [1]Sedgewick, R. 1988, Algorithms, 2nd ed. (Reading, MA: Addison-Wesley), Chapter 30.Chapter 9. Root Finding andNonlinear Sets of Equations9.0 IntroductionWe now consider that most basic of tasks, solving equations numerically. Whilemost equations are born with both a right-hand side and a left-hand side, onetraditionally moves all terms to the left, leavingf(x) = 0(9.0.1)whose solution or solutions are desired. When there is only one independent variable,the problem is one-dimensional, namely to find the root or roots of a function.With more than one independent variable, more than one equation can besatisfied simultaneously. You likely once learned the implicit function theoremwhich (in this context) gives us the hope of satisfying N equations in N unknownssimultaneously.
Note that we have only hope, not certainty. A nonlinear set ofequations may have no (real) solutions at all. Contrariwise, it may have more thanone solution. The implicit function theorem tells us that “generically” the solutionswill be distinct, pointlike, and separated from each other. If, however, life is sounkind as to present you with a nongeneric, i.e., degenerate, case, then you can geta continuous family of solutions.
In vector notation, we want to find one or moreN -dimensional solution vectors x such thatf(x) = 0(9.0.2)where f is the N -dimensional vector-valued function whose components are theindividual equations to be satisfied simultaneously.Don’t be fooled by the apparent notational similarity of equations (9.0.2) and(9.0.1).
Simultaneous solution of equations in N dimensions is much more difficultthan finding roots in the one-dimensional case. The principal difference between oneand many dimensions is that, in one dimension, it is possible to bracket or “trap” a rootbetween bracketing values, and then hunt it down like a rabbit. In multidimensions,you can never be sure that the root is there at all until you have found it.Except in linear problems, root finding invariably proceeds by iteration, andthis is equally true in one or in many dimensions.
Starting from some approximatetrial solution, a useful algorithm will improve the solution until some predeterminedconvergence criterion is satisfied. For smoothly varying functions, good algorithms347348Chapter 9.Root Finding and Nonlinear Sets of Equationswill always converge, provided that the initial guess is good enough. Indeed one caneven determine in advance the rate of convergence of most algorithms.It cannot be overemphasized, however, how crucially success depends onhaving a good first guess for the solution, especially for multidimensional problems.This crucial beginning usually depends on analysis rather than numerics. Carefullycrafted initial estimates reward you not only with reduced computational effort, butalso with understanding and increased self-esteem.
Hamming’s motto, “the purposeof computing is insight, not numbers,” is particularly apt in the area of findingroots. You should repeat this motto aloud whenever your program converges, withten-digit accuracy, to the wrong root of a problem, or whenever it fails to convergebecause there is actually no root, or because there is a root but your initial estimatewas not sufficiently close to it.“This talk of insight is all very well, but what do I actually do?” For onedimensional root finding, it is possible to give some straightforward answers: Youshould try to get some idea of what your function looks like before trying to findits roots. If you need to mass-produce roots for many different functions, then youshould at least know what some typical members of the ensemble look like. Next,you should always bracket a root, that is, know that the function changes sign in anidentified interval, before trying to converge to the root’s value.Finally (this is advice with which some daring souls might disagree, butwe give it nonetheless) never let your iteration method get outside of the bestbracketing bounds obtained at any stage.
We will see below that some pedagogicallyimportant algorithms, such as secant method or Newton-Raphson, can violate thislast constraint, and are thus not recommended unless certain fixups are implemented.Multiple roots, or very close roots, are a real problem, especially if themultiplicity is an even number. In that case, there may be no readily apparentsign change in the function, so the notion of bracketing a root — and maintainingthe bracket — becomes difficult. We are hard-liners: we nevertheless insist onbracketing a root, even if it takes the minimum-searching techniques of Chapter 10to determine whether a tantalizing dip in the function really does cross zero or not.(You can easily modify the simple golden section routine of §10.1 to return earlyif it detects a sign change in the function.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
