1*DECK POLFIT 2 SUBROUTINE POLFIT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A) 3C***BEGIN PROLOGUE POLFIT 4C***PURPOSE Fit discrete data in a least squares sense by polynomials 5C in one variable. 6C***LIBRARY SLATEC 7C***CATEGORY K1A1A2 8C***TYPE SINGLE PRECISION (POLFIT-S, DPOLFT-D) 9C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT 10C***AUTHOR Shampine, L. F., (SNLA) 11C Davenport, S. M., (SNLA) 12C Huddleston, R. E., (SNLL) 13C***DESCRIPTION 14C 15C Abstract 16C 17C Given a collection of points X(I) and a set of values Y(I) which 18C correspond to some function or measurement at each of the X(I), 19C subroutine POLFIT computes the weighted least-squares polynomial 20C fits of all degrees up to some degree either specified by the user 21C or determined by the routine. The fits thus obtained are in 22C orthogonal polynomial form. Subroutine PVALUE may then be 23C called to evaluate the fitted polynomials and any of their 24C derivatives at any point. The subroutine PCOEF may be used to 25C express the polynomial fits as powers of (X-C) for any specified 26C point C. 27C 28C The parameters for POLFIT are 29C 30C Input -- 31C N - the number of data points. The arrays X, Y and W 32C must be dimensioned at least N (N .GE. 1). 33C X - array of values of the independent variable. These 34C values may appear in any order and need not all be 35C distinct. 36C Y - array of corresponding function values. 37C W - array of positive values to be used as weights. If 38C W(1) is negative, POLFIT will set all the weights 39C to 1.0, which means unweighted least squares error 40C will be minimized. To minimize relative error, the 41C user should set the weights to: W(I) = 1.0/Y(I)**2, 42C I = 1,...,N . 43C MAXDEG - maximum degree to be allowed for polynomial fit. 44C MAXDEG may be any non-negative integer less than N. 45C Note -- MAXDEG cannot be equal to N-1 when a 46C statistical test is to be used for degree selection, 47C i.e., when input value of EPS is negative. 48C EPS - specifies the criterion to be used in determining 49C the degree of fit to be computed. 50C (1) If EPS is input negative, POLFIT chooses the 51C degree based on a statistical F test of 52C significance. One of three possible 53C significance levels will be used: .01, .05 or 54C .10. If EPS=-1.0 , the routine will 55C automatically select one of these levels based 56C on the number of data points and the maximum 57C degree to be considered. If EPS is input as 58C -.01, -.05, or -.10, a significance level of 59C .01, .05, or .10, respectively, will be used. 60C (2) If EPS is set to 0., POLFIT computes the 61C polynomials of degrees 0 through MAXDEG . 62C (3) If EPS is input positive, EPS is the RMS 63C error tolerance which must be satisfied by the 64C fitted polynomial. POLFIT will increase the 65C degree of fit until this criterion is met or 66C until the maximum degree is reached. 67C 68C Output -- 69C NDEG - degree of the highest degree fit computed. 70C EPS - RMS error of the polynomial of degree NDEG . 71C R - vector of dimension at least NDEG containing values 72C of the fit of degree NDEG at each of the X(I) . 73C Except when the statistical test is used, these 74C values are more accurate than results from subroutine 75C PVALUE normally are. 76C IERR - error flag with the following possible values. 77C 1 -- indicates normal execution, i.e., either 78C (1) the input value of EPS was negative, and the 79C computed polynomial fit of degree NDEG 80C satisfies the specified F test, or 81C (2) the input value of EPS was 0., and the fits of 82C all degrees up to MAXDEG are complete, or 83C (3) the input value of EPS was positive, and the 84C polynomial of degree NDEG satisfies the RMS 85C error requirement. 86C 2 -- invalid input parameter. At least one of the input 87C parameters has an illegal value and must be corrected 88C before POLFIT can proceed. Valid input results 89C when the following restrictions are observed 90C N .GE. 1 91C 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0. 92C 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0. 93C W(1)=-1.0 or W(I) .GT. 0., I=1,...,N . 94C 3 -- cannot satisfy the RMS error requirement with a 95C polynomial of degree no greater than MAXDEG . Best 96C fit found is of degree MAXDEG . 97C 4 -- cannot satisfy the test for significance using 98C current value of MAXDEG . Statistically, the 99C best fit found is of order NORD . (In this case, 100C NDEG will have one of the values: MAXDEG-2, 101C MAXDEG-1, or MAXDEG). Using a higher value of 102C MAXDEG may result in passing the test. 103C A - work and output array having at least 3N+3MAXDEG+3 104C locations 105C 106C Note - POLFIT calculates all fits of degrees up to and including 107C NDEG . Any or all of these fits can be evaluated or 108C expressed as powers of (X-C) using PVALUE and PCOEF 109C after just one call to POLFIT . 110C 111C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, 112C Curve fitting by polynomials in one variable, Report 113C SLA-74-0270, Sandia Laboratories, June 1974. 114C***ROUTINES CALLED PVALUE, XERMSG 115C***REVISION HISTORY (YYMMDD) 116C 740601 DATE WRITTEN 117C 890531 Changed all specific intrinsics to generic. (WRB) 118C 890531 REVISION DATE from Version 3.2 119C 891214 Prologue converted to Version 4.0 format. (BAB) 120C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 121C 920501 Reformatted the REFERENCES section. (WRB) 122C 920527 Corrected erroneous statements in DESCRIPTION. (WRB) 123C***END PROLOGUE POLFIT 124 DOUBLE PRECISION TEMD1,TEMD2 125 DIMENSION X(*), Y(*), W(*), R(*), A(*) 126 DIMENSION CO(4,3) 127 SAVE CO 128 DATA CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2), 129 1 CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3), 130 2 CO(4,3)/-13.086850,-2.4648165,-3.3846535,-1.2973162, 131 3 -3.3381146,-1.7812271,-3.2578406,-1.6589279, 132 4 -1.6282703,-1.3152745,-3.2640179,-1.9829776/ 133C***FIRST EXECUTABLE STATEMENT POLFIT 134 M = ABS(N) 135 IF (M .EQ. 0) GO TO 30 136 IF (MAXDEG .LT. 0) GO TO 30 137 A(1) = MAXDEG 138 MOP1 = MAXDEG + 1 139 IF (M .LT. MOP1) GO TO 30 140 IF (EPS .LT. 0.0 .AND. M .EQ. MOP1) GO TO 30 141 XM = M 142 ETST = EPS*EPS*XM 143 IF (W(1) .LT. 0.0) GO TO 2 144 DO 1 I = 1,M 145 IF (W(I) .LE. 0.0) GO TO 30 146 1 CONTINUE 147 GO TO 4 148 2 DO 3 I = 1,M 149 3 W(I) = 1.0 150 4 IF (EPS .GE. 0.0) GO TO 8 151C 152C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR 153C CHOOSING DEGREE OF POLYNOMIAL FIT 154C 155 IF (EPS .GT. (-.55)) GO TO 5 156 IDEGF = M - MAXDEG - 1 157 KSIG = 1 158 IF (IDEGF .LT. 10) KSIG = 2 159 IF (IDEGF .LT. 5) KSIG = 3 160 GO TO 8 161 5 KSIG = 1 162 IF (EPS .LT. (-.03)) KSIG = 2 163 IF (EPS .LT. (-.07)) KSIG = 3 164C 165C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING 166C 167 8 K1 = MAXDEG + 1 168 K2 = K1 + MAXDEG 169 K3 = K2 + MAXDEG + 2 170 K4 = K3 + M 171 K5 = K4 + M 172 DO 9 I = 2,K4 173 9 A(I) = 0.0 174 W11 = 0.0 175 IF (N .LT. 0) GO TO 11 176C 177C UNCONSTRAINED CASE 178C 179 DO 10 I = 1,M 180 K4PI = K4 + I 181 A(K4PI) = 1.0 182 10 W11 = W11 + W(I) 183 GO TO 13 184C 185C CONSTRAINED CASE 186C 187 11 DO 12 I = 1,M 188 K4PI = K4 + I 189 12 W11 = W11 + W(I)*A(K4PI)**2 190C 191C COMPUTE FIT OF DEGREE ZERO 192C 193 13 TEMD1 = 0.0D0 194 DO 14 I = 1,M 195 K4PI = K4 + I 196 TEMD1 = TEMD1 + DBLE(W(I))*DBLE(Y(I))*DBLE(A(K4PI)) 197 14 CONTINUE 198 TEMD1 = TEMD1/DBLE(W11) 199 A(K2+1) = TEMD1 200 SIGJ = 0.0 201 DO 15 I = 1,M 202 K4PI = K4 + I 203 K5PI = K5 + I 204 TEMD2 = TEMD1*DBLE(A(K4PI)) 205 R(I) = TEMD2 206 A(K5PI) = TEMD2 - DBLE(R(I)) 207 15 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 208 J = 0 209C 210C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION 211C 212 IF (EPS) 24,26,27 213C 214C INCREMENT DEGREE 215C 216 16 J = J + 1 217 JP1 = J + 1 218 K1PJ = K1 + J 219 K2PJ = K2 + J 220 SIGJM1 = SIGJ 221C 222C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1 223C 224 IF (J .GT. 1) A(K1PJ) = W11/W1 225C 226C COMPUTE NEW A COEFFICIENT 227C 228 TEMD1 = 0.0D0 229 DO 18 I = 1,M 230 K4PI = K4 + I 231 TEMD2 = A(K4PI) 232 TEMD1 = TEMD1 + DBLE(X(I))*DBLE(W(I))*TEMD2*TEMD2 233 18 CONTINUE 234 A(JP1) = TEMD1/DBLE(W11) 235C 236C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS 237C 238 W1 = W11 239 W11 = 0.0 240 DO 19 I = 1,M 241 K3PI = K3 + I 242 K4PI = K4 + I 243 TEMP = A(K3PI) 244 A(K3PI) = A(K4PI) 245 A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP 246 19 W11 = W11 + W(I)*A(K4PI)**2 247C 248C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE 249C PRECISION 250C 251 TEMD1 = 0.0D0 252 DO 20 I = 1,M 253 K4PI = K4 + I 254 K5PI = K5 + I 255 TEMD2 = DBLE(W(I))*DBLE((Y(I)-R(I))-A(K5PI))*DBLE(A(K4PI)) 256 20 TEMD1 = TEMD1 + TEMD2 257 TEMD1 = TEMD1/DBLE(W11) 258 A(K2PJ+1) = TEMD1 259C 260C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND 261C ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE 262C COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT, 263C THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST 264C SIGNIFICANT BITS ARE IN A(K5PI) . 265C 266 SIGJ = 0.0 267 DO 21 I = 1,M 268 K4PI = K4 + I 269 K5PI = K5 + I 270 TEMD2 = DBLE(R(I)) + DBLE(A(K5PI)) + TEMD1*DBLE(A(K4PI)) 271 R(I) = TEMD2 272 A(K5PI) = TEMD2 - DBLE(R(I)) 273 21 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 274C 275C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE 276C MAXDEG HAS BEEN REACHED 277C 278 IF (EPS) 23,26,27 279C 280C COMPUTE F STATISTICS (INPUT EPS .LT. 0.) 281C 282 23 IF (SIGJ .EQ. 0.0) GO TO 29 283 DEGF = M - J - 1 284 DEN = (CO(4,KSIG)*DEGF + 1.0)*DEGF 285 FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN 286 FCRIT = FCRIT*FCRIT 287 F = (SIGJM1 - SIGJ)*DEGF/SIGJ 288 IF (F .LT. FCRIT) GO TO 25 289C 290C POLYNOMIAL OF DEGREE J SATISFIES F TEST 291C 292 24 SIGPAS = SIGJ 293 JPAS = J 294 NFAIL = 0 295 IF (MAXDEG .EQ. J) GO TO 32 296 GO TO 16 297C 298C POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE 299C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND. 300C 301 25 NFAIL = NFAIL + 1 302 IF (NFAIL .GE. 3) GO TO 29 303 IF (MAXDEG .EQ. J) GO TO 32 304 GO TO 16 305C 306C RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT 307C EPS = 0.) 308C 309 26 IF (MAXDEG .EQ. J) GO TO 28 310 GO TO 16 311C 312C SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.) 313C 314 27 IF (SIGJ .LE. ETST) GO TO 28 315 IF (MAXDEG .EQ. J) GO TO 31 316 GO TO 16 317C 318C RETURNS 319C 320 28 IERR = 1 321 NDEG = J 322 SIG = SIGJ 323 GO TO 33 324 29 IERR = 1 325 NDEG = JPAS 326 SIG = SIGPAS 327 GO TO 33 328 30 IERR = 2 329 CALL XERMSG ('SLATEC', 'POLFIT', 'INVALID INPUT PARAMETER.', 2, 330 + 1) 331 GO TO 37 332 31 IERR = 3 333 NDEG = MAXDEG 334 SIG = SIGJ 335 GO TO 33 336 32 IERR = 4 337 NDEG = JPAS 338 SIG = SIGPAS 339C 340 33 A(K3) = NDEG 341C 342C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT 343C ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES 344C 345 IF(EPS .GE. 0.0 .OR. NDEG .EQ. MAXDEG) GO TO 36 346 NDER = 0 347 DO 35 I = 1,M 348 CALL PVALUE (NDEG,NDER,X(I),R(I),YP,A) 349 35 CONTINUE 350 36 EPS = SQRT(SIG/XM) 351 37 RETURN 352 END 353