Nash - Compact Numerical Methods for Computers (523163), страница 30
Текст из файла (страница 30)
Cautious users may wish to pass back thenumber of sweeps and the limit value and perform a test here.}{STEP 3 is not needed in this version of the algorithm.}{STEP 4 -- transformation of initial eigenvectors and creation of matrixC, which we store in array B}for i := 1 to n dobeginif B[i,i]c=0.0 then halt; {matrix B is not computationallypositive definite.}s := 1.0/sqrt(B[i,i]);for j := 1 to n do V[j,i] := s * V[j,i]; {to form Bihalf}end; {loop on i}{STEP 5 -- not needed as matrix A already entered}{STEP 6 -- perform the transformation 11.11b}for i := l to n dobeginfor j := i to n do {NOTE: i to n NOT 1 to n}begins := 0.0;for k := 1 to n dofor m := 1 to n dos := s+V[k,i]*A[k,m]*V[m,j];B[i.j] := s; B[j,i] := s;end; {loop on j}end; {loop on i}{STEP 7 -- revised to provide simpler Pascal code}initev := false; {Do not initialize eigenvectors, since we haveprovided the initialization in STEP 4.}writeln(‘Eigensolutions of general problem’);evJacobi(n, B, V, initev); {eigensolutions of generalized problem}end; {alg15.pas == genevJac}Example 11.1.
The generalised symmetric eigenproblem: the anharmonic oscillatorThe anharmonic oscillator problem in quantum mechanics provides a system forwhich the stationary states cannot be found exactly. However, they can becomputed quite accurately by the Rayleigh-Ritz method outlined in example 2.5.The Hamiltonian operator (Newing and Cunningham 1967) for the anharmonicoscillator is writtenH = –d2/dx 2+ k 2 x 2+ k 4 x 4 .(11.22)The eigenproblem arises by minimising the Rayleigh quotient(11.23)where φ(x ) is expressed as some linear combination of basis functions.
If thesebasis functions are orthonormal under the integration over the real line, then aconventional eigenproblem results; However, it is common that these functions139The generalised symmetric matrix eigenvalue problemare not orthonormal, so that a generalised eigenproblem arises. Because of thenature of the operations these eigenproblems are symmetric.In order to find suitable functions fj(x) to expand φ as(11.24)we note that the oscillator which has coefficients k 2 = 1, k 4= 0 in its potential hasexact solutions which are polynomials multiplied byexp(-0·5 x 2 ) .Therefore, the basis functions which will be used here aref i (x) = Nxj -1 exp(-α x 2 )(11.25)where N is a normalising constant.The approximation sought will be limited to n terms.
Note now thatHfj (x) = N exp(–α x 2)[–(j – 1)(j – 2)x j-3 + 2α (2j – 1)x j- 1+( k 2 -4a2 )x j +1+ k4 xj+3].(11.26)The minimisation of the Rayleigh quotient with respect to the coefficients cj givesthe eigenproblem(11.27)Ac = eBcwhere(11.28)and(11.29)These integrals can be decomposed to give expressions which involve only theintegralsfor m oddfor m even(11.30)for m = 0.The normalising constant N 2 has been chosen to cancel some awkard constants inthe integrals (see, for instance, Pierce and Foster 1956, p 68).Because of the properties of the integrals (11.30) the eigenvalue problem(11.27) reduces to two smaller ones for the even and the odd functions.
If we set aparity indicator w equal to zero for the even case and one for the odd case,140Compact numerical methods for computerswe can substitutej - 1 = 2 (q -1) + w(11.31 a)(11.31b)i- 1 = 2 (p-1) + wwhere p and q will be the new indices for the matrices A and B running from 1 ton'= n /2 (assuming n even). Thus the matrix elements areÃp q=-(j – 1)(j – 2)Is + 2 a(2j – 1)I s+2 + (k 2– 4a 2 )Is +4+ k4 I s + 6(11.32)and(11.33)wheres = i + j – 4 =2(p + q – 3+ w)and j is given by (11.31a). The tilde is used to indicate the re-numeration of Aand B.The integrals (11.30) are easily computed recursively.STEP0123DESCRIPTIONEnter s, α. Note s is even.Let v = 1.If s<0, stop.
Is is in v. For s<0 this is always multiplied by 0.For k = 1 to s/2.Let v = v * (2 * k- 1) * 0·25/ α.End loop on k.4End integral. Is is returned in v.As an example, consider the exactly solvable problem using n' = 2, and α = 0·5for w = 0 (even parity). Then the eigenproblem haswith solutionse =1c = (1, 0)Te=5c = 2-½(-1, 2)T.andThe same oscillator (a = 0·5) with w = 1 and n' = 10 should also have exactsolutions. However, the matrix elements range from 0·5 to 3·2E+17 and thesolutions are almost all poor approximations when found by algorithm 15.Likewise, while the problem defined by n' = 5, w = 0, a = 2, k2= 0, k4 = 1 issolved quite easily to give the smallest eigenvalue e1= 1·06051 with eigenvectorc = (0·747087, 1·07358, 0·866449, 0·086206, 0·195257)Tthe similar problem with n' = 10 proves to have a B matrix which is computationally singular (step 4 of algorithm 15).
Inverse iteration saves the day giving, forn ' = 5, e = 1·0651 and for n' = 10, e = 1·06027 with eigenvectors having smallresiduals. These results were found using an inverse iteration program based onThe generalised symmetric matrix eigenvalue problem141Gauss elimination run on a Data General NOVA in 23-bit binary arithmetic.The particular basis functions (11.25) are, as it turns out, rather a poor set tochoose since no matter what exponential parameter α is chosen, the matrix B asgiven by (11.29) suffers large differences in scale in its elements. To preserve thesymmetry of the problem, a diagonal scaling could be applied, but a better idea isto use basic functions which are already orthogonal under the inner productdefined by (11.29).
These are the Hermite functions (Newing and Cunningham1967). As the details of the solution are not germane to the discussion, they willnot be given.BASIC program code to compute the matrix elements for this problem has been givenin Nash (1984a, §18-4, in particular the code on p 242). The reader is warned thatthere are numerous misprints in the referenced work, which was published about thetime of the failure of its publisher. However, the BASIC code, reproduced photographically, is believed to be correct.Previous HomeChapter 12OPTIMISATION AND NONLINEAR EQUATIONS12.1. FORMAL PROBLEMS IN UNCONSTRAINEDOPTIMISATION AND NONLINEAR EQUATIONSThe material which follows in the next few chapters deals with finding the minimaof functions or the roots of equations. The functions or equations will in generalbe nonlinear in the parameters, that is following Kowalik and Osborne (1968),the problems they generate will not be solvable by means of linear equations(though, as we shall see, iterative methods based on linear subproblems are veryimportant).A special case of the minimisation problem is the nonlinear least-squaresproblem which, because of its practical importance, will be presented first.
Thiscan be stated: given M nonlinear functions of n parametersi = 1, 2, . . . , Mf i(b 1 ,b2, . . . ,b n)(12.1)minimise the sum of squares(12.2)It is convenient to collect the n parameters into a vector b; likewise the functionscan be collected as the vector of M elements f. The nonlinear least-squaresproblem commonly, but not exclusively, arises in the fitting of equations to databy the adjustment of parameters.
The data, in the form of K variables (whereK = 0 when there are no data), may be thought of as occupying a matrix Y ofwhich the j th column yj gives the value of variable j at each of the M data points.The terms parameters, variables and data points as used here should be noted,since they lead naturally to the least-squares problem with(12.3)f( b , Y) = g(b, yl, y2, . . . , yK -1) – yKin which it is hoped to fit the function(s) g to the variable yK by adjusting theparameters b. This will reduce the size of the functions f, and hence reduce thesum of squares S.
Since the functions f in equation (12.3) are formed asdifferences, they will be termed residuals. Note that the sign of these is theopposite of that usually chosen. The sum of squares is the same of course. Mypreference for the form used here stems from the fact that the partial derivativesof f with respect to the parameters b are identical to those of g with respect to thesame parameters and the possibility of the sign error in their computation isavoided.142NextOptimisation and nonlinear equations143By using the shorthand of vector notation, the nonlinear least-squares problemis written: minimiseS ( b) = f Tf = f T (b, Y)f(b, Y)(12.4)with respect to the parameters b.
Once again, K is the number of variables, M isthe number of data points and n is the number of parameters.Every nonlinear least-squares problem as defined above is an unconstrainedminimisation problem, though the converse is not true. In later sections methodswill be presented with aim to minimise S(b) where S is any function of theparameters. Some ways of handling constraints will also be mentioned. Unfortunately, the mathematical programming problem, in which a minimum is sought fora function subject to many constraints, will only be touched upon briefly sincevery little progress has been made to date in developing algorithms with minimalstorage requirements.There is also a close relationship between the nonlinear least-squares problemand the problem of finding solutions of systems of nonlinear equations.
A systemof nonlinear equationsf (b,Y) = 0(12.5)having n = M (number of parameters equal to the number of equations) can beapproached as the nonlinear least-squares problem: minimiseS = f Tf(12.6)with respect to b. For M greater than n, solutions can be sought in the leastsquares sense; from this viewpoint the problems are then indistinguishable. Theminimum in (12.6) should be found with S = 0 if the system of equations has asolution. Conversely, the derivativesj = 1, 2, . . . , n(12.7)for an unconstrained minimisation problem, and in particular a least-squaresproblem, should be zero at the minimum of the function S( b), so that theseproblems may be solved by a method for nonlinear equations, though localmaxima and saddle points of the function will also have zero derivatives and areacceptable solutions of the nonlinear equations. In fact, very little research hasbeen done on the general minimisation or nonlinear-equation problem whereeither all solutions or extrema are required or a global minimum is to be found.The minimisation problem when n = 1 is of particular interest as a subproblemin some of the methods to be discussed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
