Nash - Compact Numerical Methods for Computers (523163), страница 18
Текст из файла (страница 18)
. . , n, is in the diagonal or pivot position. Thismodified procedure, called Gauss elimination with partial pivoting, has a largeliterature (see Wilkinson (1965) for a discussion with error analysis). Since therows of A have been exchanged, this procedure gives the triangular decompositionof a transformed matrixPA = LR(6.27)where P is the permutation matrix, simply a re-ordering of the columns of a unitmatrix appropriate to the re-ordering of the rows of A.Particular methods exist for further minimising error propagation in the Gausselimination procedure by using double-precision accumulation of vector innerproducts. These go by the names of Crout and Doolittle and are discussed, forinstance, by Dahlquist and Björck (1974) as well as by Wilkinson (1965).
Since theaccumulation of inner products in several computer programming languages on avariety of machines is a non-trivial operation (though on a few machines it issimpler to accumulate in double than in single precision), these will not bediscussed here. Nor will the Gauss elimination with complete pivoting, whichchooses as pivot the largest element in the current matrix, thereby requiring bothcolumn and row permutations. The consensus of opinion in the literature appearsto be that this involves too much work for the very few occasions when completepivoting is distinctly more accurate than partial pivoting.The Gauss elimination with partial pivoting and back-substitution are nowstated explicitly.
All work is done in an array A of dimension n by n + p, where pis the number of right-hand sides b to be solved.Algorithm 5. Gauss elimination with partial pivotingProcedure gelim ( n : integer; {order of equations}p : integer; {number of right hand sides}var A : rmatrix; {equation coefficients in row order with righthand sides built into the matrix columns n+1, . . ,n+p}to1 : real); {pivot tolerance}{alg05.pas == Gauss elimination with partial pivoting.This form does not save the pivot ordering, nor does it keepthe matrix decomposition which results from the calculations.Copyright 1988 J.
C. Nash}vardet, s : real;h,i,j,k: integer;begin {STEP 0}det := 1.0; {to initialise determinant value}writeln(‘alg05.pas -- Gauss elimination with partial pivoting’);for j := 1 to (n-1) do {STEP 1}76Compact numerical methods for computersAlgorithm 5. Gauss elimination with partial pivoting (cont.)begin {STEP 2a}s := abs(A[j,j]); k := j;for h := (j+l) to n do {STEP 2b}beginif abs(A[h,j])>s thenbegins := abs(A[h,j]); k := h;end;end; {loop on h}if k<>j then {STEP 3 -- perform row interchange. Here we do thisexplicitly and trade time for the space to store indicesand the more complicated program code.}beginwriteln(‘Interchanging rows ’,k,‘ and ’,j);for i := j to (n+p) dobegins := A[k,i]; A[k,i] := A[j,i]; A[j,i] := s;end; {loop on i}det := det; {to change sign on determinant because of interchange}end; (interchange)det := det*A[j,j]; {STEP 4}if abs(A[j,j])<tol thenbeginwriteln(‘Matrix computationally singular -- pivot < ’,tol);halt;end;for k := (j+l) to n do {STEPS}beginA[k,j] := A[k,j]/A[j,j]; {to form multiplier m[k,j]}for i := (j+1) to (n+p) doA[k,i] := A[k,i]-A[k,j]*A[j,i]; {main elimination step}end; {loop on k -- STEP 6}det := det*A[n,n]; {STEP 7}if abs(A[n,n])<tol thenbeginwriteln(‘Matrix computationally singular -- pivot < ’,tol);halt;end;end, {loop on j -- this ends the elimination steps}writeln(‘Gauss elimination complete -- determinant = ’,det);end; {alg05.pas}The array A now contains the informationmij = A[i, j]for i > jfor i < j < nRij = A[i, j]fi(j) = A[i,j]for j = n+1, n+2, .
. . , n+pthat is, the elements of the p right-hand sides.Linear equations—a direct approach77Algorithm 6. Gauss elimination back-substitutionThis algorithm is designed to follow Gauss elimination (algorithm 5), but can be applied alsoto systems of equations which are already triangular or which have been brought to triangularform by other methods, for instance, the Givens’ reduction of §4.2 (algorithm 3).procedure gebacksub(n, p:integer; {size of problem n=nRow, p=nRHS}var A : rmatrix); {work array containingGauss elimination reduced coefficientmatrix and transformed right hand sides}{alg06.pas == Gauss elimination back-substitution.Places solutions to linear equation systems in columns n+1,..,n+pof matrix A.
Alg05.pas (Gauss elimination) must be executed first withright hand sides in the columns n+1,..,n+p of matrix A in orderto triangularize the system of equations.Copyright 1988 J. C. Nash}vars : real; {accumulator}i, j, k: integer;beginwriteln(‘alg06.pas -- Gauss elimination back-substitution’);for i:=(n+1) to (n+p) do {STEP 1}beginA[n,i]:=A[n,i]/A[n,n]; {STEP 2}for j:=(n-1) down to 1 do {STEP 3}begins:=A[j,i]; {STEP 4}for k:=(j+1) to n do {STEP 5}begins:=s-A[j,k]*A[k,i]; {to subtract contributions from solutionelements which have already been determined}end; {loop on k}A[j,i]:=s/A[j,j]; {STEP 6 -- to fix solution element j}end; {loop on j -- STEP 7}end; {loop on i -- STEP 8}end; {alg06.pas}The solutions to the triangular system(s) of equations and hence to the original equations(2.2) are contained in columns n + 1, n + 2, . .
. , n+p, of the working array.Example 6.1. The use of linear equations and linear least-squares problemsOrganisations which publish statistics frequently use indices to summarise thechange in some set of measurable quantities. Already in example 3.2 we haveused indices of the use of various chemicals in agriculture and an index for farmincome. The consumer price index, and the Dow Jones and Financial Timesindices provide other examples.
Such indices are computed by dividing theaverage value of the quantity for period t by the average for some base periodt = 0 which is usually given the index value 100. Thus, if the quantity is called P,then(6.28)78Compact numerical methods for computerswhere(6.29)given n classes or types of quantity P, of which the jth has value Ptj in period tand is assigned weight Wj in the average.
Note that it is assumed that theweighting Wj is independent of the period, that is, of time. However, theweightings or ‘shopping basket’ may in fact change from time to time to reflectchanging patterns of product composition, industrial processes causing pollution,stocks or securities in a portfolio, or consumer spending.Substitution of (6.29) into (6.28) givesFinally, letting(6.30)gives(6.31)Thus, if n periods of data It, Ptj, j = 1, . .
. , n, are available, we can compute theweightings KWj. Hence, by assuming(6.32)that is, that the weights are fractional contributions of each component, we canfind the value of K and each of the Wj. This involves no more nor less than thesolution of a set of linear equations. The work of solving these is, of course,unnecessary if the person who computes the index publishes his set of weights-asindeed is the case for several indices published in the Monthly Digest of Staristics † .Unfortunately, many workers do not deem this a useful or courteous practicetowards their colleagues, and I have on two occasions had to attempt to discoverthe weightings.
In both cases it was not possible to find a consistent set of weightsover more than n periods, indicating that these were being adjusted over time.This created some difficulties for my colleagues who brought me the problems,since they were being asked to use current price data to generate a provisionalestimate of a price index considerably in advance of the publication of the indicesby the agency which normally performed the task. Without the weights, or evenapproximate values from the latest period for which they were available, it wasnot possible to construct such estimates.
In one case the calculation was to have† Monthly Digest of Statistics UK Central Statistical Office (London: HMSO).Linear equations—a direct approach79TABLE 6.1. Prices and indices.P1P11·11·ll·25l·31·281·310·50·50·50·60·60·60·62PP 43l·31·36l·41·411·4121·521·63·63·63·63·63·953·93·95I1100103·718104·487109·167114·974115·897118·846I2100103·718104·487109·167114·97498·4615101·506used various proposed oil price levels to ascertain an index of agricultural costs.When it proved impossible to construct a set of consistent weights, it wasnecessary to try to track down the author of the earlier index values.As an example of such calculations, consider the set of prices shown in table 6.1and two indices I1 and I2 calculated from them. I1 is computed using proportions0·4, 0·1, 0·3 and 0·2 respectively of P1, P2, P3 and P4.
I2 uses the same weightsexcept for the last two periods where the values 0·35, 0·15, 0·4 and 0·1 are used.Suppose now that these weights are unknown. Then the data for the first fourperiods give a set of four equations (6.31) which can be solved to giveKW =using Gauss elimination (Data General NOVA, 23-bit binary mantissa). Applyingthe normalisation (6.32) givesW =If these weights are used to generate index numbers for the last three periods, thevalues I1 will be essentially reproduced, and we would detect a change in theweighting pattern if the values I2 were expected.An alternative method is to use a least-squares formulation, since if the set ofweights is consistent, the residual sum of squares will be zero.
Note that there isno constant term (column of ones) in the equations. Again on the NOVA in23-bit arithmetic, I1 gives80Compact numerical methods for computerswith a residual sum of squares (using KW) over the seven periods of 4·15777E–7.The same calculation with I2 gives a residual sum of squares of 241·112, showingthat there is not a consistent set of weights.
It is, of course, possible to find aconsistent set of weights even though index numbers have been computed using avarying set; for instance, if our price data had two elements identical in oneperiod, any pair of weights for these prices whose sum was fixed would generatethe same index number.6.3. VARIATIONS ON THE THEME OF GAUSS ELIMINATIONGauss elimination really presents only one type of difficulty to the programmerwhich of the many possible variations to implement. We have already touchedupon the existence of two of the better known ones, those of Crout and Doolittle(see Dahlquist and Björck 1974, pp 157–8). While these methods are useful andimportant, they require double-precision arithmetic to be used to full advantage,so cannot be used effectively if the computing system at hand lacks this capability.Bowdler et al (1966) present ALGOL versions of the Crout algorithm whichimplicitly scale the rows of the coefficient matrix.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
