Abstract

            This document contains a fit to AiGauss(z) = erfc(z)/2 with an error ε < 10-15 × AiGauss(z).  #BestFit

The code is in ..\..\zips\FittingBrack.zip#ZipFile is an annotated directory of this file.

AiGauss

The function is defined by

(1.1)

This is approximated as

(1.2)

Where

(1.3)

The form in (1.2) has 15 digit accuracy for the exponential and also for the Pade approximate.  This is better than the form Pade = log( AiGau(z)) which has 15 digit accuracy in the log, but uses two of these digits at z = -6 where the lot is approximately -36 to simply represent the order of magnitude of the function.

 

The first fit showed errors in the integral version of AiGauss itself ../aigau15/Fitting%20Brack.doc .htm

 

The code for producing AiGauss data is saved here as ..\..\zips\genint.zip  Compiling was done with

Wfl386 tgleg67

This produced a welcome extra digit of accuracy beyond that available with Watfor. 

 

The end

Aigau2.dir

AIGAU.DIR, OUTPUT FILE NAME

..\TestData\BRACK.OUT,    DATA SET NAME

1,      DIMENSION OF X IN F(X)

8,8 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS

0,  0 fits to exp(-x**2) times a Pade, 0 to a Pade

1E-25,1E-15,  err(i)= max(ep1,ep2*dat(i))   this fit is to a part in 100000,      NUMBER OF MINIMIZATION STEPS   1015  real*8 limit.

 

.9, INITIAL FR DECREASE DESIRED   as we lowe precision small jumps are harder than big ones

1 VARY O FIX   FITTED CONS  ERROR IN CONS

    0  0.5000000000000000       +- 0.000E+00

    1 -0.7780697362924373       +- 0.916E-04

    1  0.5990638419338665       +- 0.129E-03

    1 -0.2771270371606460       +- 0.891E-04

    1  0.8027282341142808E-01   +- 0.354E-04

    1 -0.1379770398066123E-01   +- 0.805E-05

    1  0.1113061469828807E-02   +- 0.857E-06

    1 0

    1  -2.684518639716602       +- 0.183E-03

    1   3.227282589432518       +- 0.465E-03

    1  -2.263586666841055       +- 0.520E-03

    1   1.006883612355409       +- 0.330E-03

    1 -0.2865296587487073       +- 0.127E-03

    1  0.4891172945745505E-01   +- 0.285E-04

    1 -0.3945696943660320E-02   +- 0.304E-05

    1 0

 CHIS=0.3550E+08 FOR  5001 DATA POINTS

 CHI USED IN CALCULATING ERRORS IS 0.3550E+08

Figure 1 fi-fA(xi)

            In this case fi is exp(z2) × AiGauss(z), and fA(x) = Brack(x) with the 8 numerator and denominator constants shown above.  The quantitiy minimized emphasizes the left side of figure 1.  The numerical thickness in the data to the left indicates how far the fit has to go to get to the noise level.

 

I:\PUBLIC~1\nlfit\EXPTIM~1>mpnlfit aigau2.dir  the complete output is below, but this is

 NADIR=aigau2.dir                               a good time for coffee.

 NAOUT=AIGAU.DIR

 NADAT=..\TestData\BRACK.OUT

 ND=           1

 VORF   CONS

    0  0.500000

    1 -0.778070

    1  0.599064

    1 -0.277127

    1  0.802728E-01

    1 -0.137977E-01

    1  0.111306E-02

    1  0.000000E+00

    1  -2.68452

    1   3.22728

    1  -2.26359

    1   1.00688

    1 -0.286530

    1  0.489117E-01

    1 -0.394570E-02

    1  0.000000E+00

 FIT IS BEING MADE TO DATA IN FILE

 ..\TestData\BRACK.OUT

FR,CHI,CHB,CHL 0.900000000      35499382.3     0.100000000E+67 0.500000000E+67

FR,CHI,CHB,CHL 0.729000000      25968696.2      25879049.7      35499382.3

FR,CHI,CHB,CHL 0.387420489      10184521.6      10060805.0      25968696.2

FR,CHI,CHB,CHL 0.581497370E-01  556637.656      592227.255      10184521.6

 from bnraph going low   cannot get to .01 * chi with any lambda  exp(-175)

FR,CHI,CHB,CHL 0.100000000E-01  16867.8334      6673.35984      556637.656

 from bnraph going low

FR,CHI,CHB,CHL 0.100000000E-01  9421.78133      9466.63002      16867.8334

 from bnraph going low

FR,CHI,CHB,CHL 0.100000000E-01  9411.07375      9407.58683      9421.78133

 from bnraph going low

FR,CHI,CHB,CHL 0.109000000      9423.80202      9409.60583      9411.07375

IREV=    1 IEND=99993

 from bnraph going low 

FR,CHI,CHB,CHL 0.910900000      9423.80202      9409.60583      9411.07375

IREV=    2 IEND=99993

 from bnraph going low

FR,CHI,CHB,CHL 0.991090000      9423.80202      9409.60583      9411.07375

IREV=   20 IEND=99993   I modified the code, so that if it keeps getting the same value                                   it jumps to the last IREV  within 20-100, the chi-square                                at this level is dominated by truncation errors.

AT END FR, CHI  0.999109000      9411.07375

 TEMP2=  0.000000000E+00

 NEW SECTION FR=        0.9991090000000000   nlfit is set to automatically restart

FR,CHI,CHB,CHL 0.999999990      9411.07375     0.100000000E+67 0.500000000E+67

FR,CHI,CHB,CHL 0.999999970      9411.07375      9411.07347      9411.07375

IREV=    1 IEND=99992

FR,CHI,CHB,CHL 0.999999997      9411.07375      9411.07372      9411.07375

IREV=    2 IEND=99992

FR,CHI,CHB,CHL 0.999999999      9411.07375      9411.07374      9411.07375

 AT END FR, CHI  0.999999999      9411.07375

 TEMP2=  0.999999901E-09

 NEW SECTION FR=        0.9999999990000000

 CHI AT END OF FIT    9411.0737524199120000

           0       0.5000000000000000 +-       0.0000000000000000

           1      -0.8591435477130216 +-  1.6272308649864600D-004

           1       0.7341466366171062 +-  2.5735696428003250D-004

           1      -0.3871217792735996 +-  2.0090444822930290D-004

           1       0.1342351250044111 +-  9.4349896125745370D-005

           1      -0.0304744298576865 +-  2.7806336792956480D-005

           1       0.0042054633161981 +-  4.8818201328104450D-006

           1 -2.7497424824452200D-004 +-  4.0523682417858010D-007

           1      -2.8466662625211940 +-  3.2544617298786200D-004

           1       3.6804121795519390 +-       0.0008819329591858

           1      -2.8327305035912660 +-       0.0010715148057411

           1       1.4258647627189870 +-       0.0007606586755338

           1      -0.4833027070743471 +-  3.4312244331430270D-004

           1       0.1085165756053468 +-  9.9288237674976140D-005

           1      -0.0149079729638434 +-  1.7305634415493940D-005

           1       0.0009747584530339 +-  1.4365264327455150D-006

 FINAL CHIT=    9411.0737524199120000

 CHI/DEG OF FREE        1.8874997497833760

 RESULTS IN FILE AIGAU.DIR

 

The entire run takes about  hour on a 1 Giga-Hertz laptop.

Aigau.dir

aigau2.dir, OUTPUT FILE NAME

..\TestData\BRACK.OUT,    DATA SET NAME                        

1,      DIMENSION OF X IN F(X)                                 

8,8 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS                

0,  0 fits to exp(-x**2) times a Pade, 0 to a Pade             

1E-25,1E-15,  err(i)= max(ep1,ep2*dat(i))                      

100000,      NUMBER OF MINIMIZATION STEPS                      

.9, INITIAL FR DECREASE DESIRED                                 

1 VARY O FIX   FITTED CONS  ERROR IN CONS                      

    0  0.5000000000000000       +- 0.000E+00

    1 -0.8591435477130216       +- 0.163E-03

    1  0.7341466366171062       +- 0.257E-03

    1 -0.3871217792735996       +- 0.201E-03

    1  0.1342351250044111       +- 0.943E-04

    1 -0.3047442985768653E-01   +- 0.278E-04

    1  0.4205463316198079E-02   +- 0.488E-05

    1 -0.2749742482445220E-03   +- 0.405E-06

    1  -2.846666262521194       +- 0.325E-03

    1   3.680412179551939       +- 0.882E-03

    1  -2.832730503591266       +- 0.107E-02

    1   1.425864762718987       +- 0.761E-03

    1 -0.4833027070743471       +- 0.343E-03

    1  0.1085165756053468       +- 0.993E-04

    1 -0.1490797296384342E-01   +- 0.173E-04

    1  0.9747584530339248E-03   +- 0.144E-05

 CHIS= 9411.     FOR  5001 DATA POINTS

 CHI USED IN CALCULATING ERRORS IS  9411.   

Figure 2  (f-fA)/(1e-15*f)

            The region -6<z<4 is approaching its accuracy limit.  This is probably caused by truncation errors in evaluating the integrals.  To the extent that these effects are “random” the underlying curve is more accurate by a factor of (Nd/Ncons) = 17.  Both regions have systematic errors on the order of 2 implying that this fit does not quite reach the goal.

            The pass with 9,9 unlike the earlier ones, does not drop chi-square by orders of magnitude.  The numerical noise ~ 100 also stops a systematic approach to the minimum, once its vicinity has been reached.

FIT IS BEING MADE TO DATA IN FILE

 ..\TestData\BRACK.OUT

FR,CHI,CHB,CHL 0.900000000      9612.87997     0.100000000E+67 0.500000000E+67

FR,CHI,CHB,CHL 0.729000000      7079.62515      7007.78950      9612.87997

 from bnraph going low

FR,CHI,CHB,CHL 0.387420489      6713.67300      6599.35985      7079.62515

 from bnraph going low

FR,CHI,CHB,CHL 0.152334763      6711.56567      6691.10203      6713.67300

 from bnraph going low

FR,CHI,CHB,CHL 0.120885292      6689.95567      6705.41629      6711.56567

 from bnraph going low

FR,CHI,CHB,CHL 0.100000000E-01  6677.24298      6685.54279      6689.95567

 from bnraph going low

FR,CHI,CHB,CHL 0.100000000E-01  6639.01696      6674.19657      6677.24298

 from bnraph going low

FR,CHI,CHB,CHL 0.100000000E-01  6694.62463      6635.20438      6639.01696

IREV=    1 IEND=99993

 from bnraph going low

FR,CHI,CHB,CHL 0.901000000      6694.62463      6635.20438      6639.01696

IREV=    2 IEND=99993

 from bnraph going low

FR,CHI,CHB,CHL 0.990100000      6694.62463      6635.20438      6639.01696

IREV=   20 IEND=99993

FR,CHI,CHB,CHL 0.999010000      6639.01696      6639.01696      6639.01696

 AT END FR, CHI  0.999010000      6639.01696

 TEMP2=  0.000000000E+00

 NEW SECTION FR=        0.9990100000000001

FR,CHI,CHB,CHL 0.999999990      6639.01696     0.100000000E+67 0.500000000E+67

FR,CHI,CHB,CHL 0.999999970      6639.01696      6639.01676      6639.01696   At this IREV=    1 IEND=99992               level there are no small predictions  rather than

FR,CHI,CHB,CHL 0.999999997      6639.01696      6639.01694      6639.01696  the

IREV=    2 IEND=99992               traditional fr  1, would do better with fr  .9 so

FR,CHI,CHB,CHL 0.999999999      6639.01696      6639.01695      6639.01696 that a larger

 AT END FR, CHI  0.999999999      6639.01696             distance in parameter space is

 TEMP2=  0.999999953E-09                                 attempted.

 NEW SECTION FR=        0.9999999990000000

Figure 3 (fi-fA(xi))/(fi × 1e-15)

            The systematic errors appear to be on the order of 5 × 10-16.  The “random” errors result in function errors that are reduced from the fluctuationa above by a factor approximately equal to 5000/18 = 16 so that these are also less than 1 × 10-15

             Following is last constant file with all updates     

aigau.dir, OUTPUT FILE NAME

..\TestData\brack.OUT,    DATA SET NAME                        

1,      DIMENSION OF X IN F(X)                                 

9,9 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS               

0,  0 fits to exp(-x**2) times a Pade, 0 to a Pade             

1E-30,1E-15,  err(i)= max(ep1,ep2*dat(i))                      

100000,      NUMBER OF MINIMIZATION STEPS                      

.9, INITIAL FR DECREASE DESIRED                                

1 VARY O FIX   FITTED CONS  ERROR IN CONS                      

    0  0.5000000000000000       +- 0.000E+00

    1  -1.065781843119121       +- 0.254E-01

    1   1.098849939447941       +- 0.441E-01

    1 -0.7059835832256011       +- 0.380E-01

    1  0.3064328550407249       +- 0.202E-01

    1 -0.9175866513292681E-01   +- 0.709E-02

    1  0.1853701251033871E-01   +- 0.163E-02

    1 -0.2323566331513835E-02   +- 0.228E-03

    1  0.1400748062111785E-03   +- 0.152E-04

    1  -3.259942853333778       +- 0.507E-01

    1   4.876151480518089       +- 0.145   

    1  -4.406424837428250       +- 0.189   

    1   2.661133284964304       +- 0.147   

    1  -1.118888036108821       +- 0.746E-01

    1  0.3293940293195529       +- 0.256E-01

    1 -0.6596029956005756E-01   +- 0.581E-02

    1  0.8236827162542478E-02   +- 0.810E-03

    1 -0.4965522790119376E-03   +- 0.540E-04

 CHIS= 6368.     FOR  5001 DATA POINTS

 CHI USED IN CALCULATING ERRORS IS  6368.   

 

            These constants above give

  

Note that the signs of the constants are such that for negative z, every term adds.  There is no cancellation and the truncation error is minimized.  There is a slight gain in fitting the results of the integral directly.  The constants below are slightly better.

Exp(-z2)Brack(z)

            The code was designed to fit AiGauss directly.  Starting with the constants above, a number of formats needed changing in ..\TestData\genint\TGLEG67.FOR to get the full range of the data  The data above with much smaller derivatives seems unaffected by these changes.  The final fit begins with the constants from fitting brack above and a chi-square of 4647.  This is lower than that for brack.  The values used now are direct integration values while brack’s fit multiplied these values by exp(z2).  It appears that this introduces an extra 2000 into the chi-square.

aigau2.dir, OUTPUT FILE NAME

..\TestData\aigauss.OUT,    DATA SET NAME                      

1,      DIMENSION OF X IN F(X)                                 

9,9 NUMBER OF NUMERATOR TERMS, DENOMINATOR TERMS               

1,  1 fits to exp(-x**2) times a Pade, 0 to a Pade             

1E-30,1E-15,  err(i)= max(ep1,ep2*dat(i))                      

100000,      NUMBER OF MINIMIZATION STEPS                      

.9, INITIAL FR DECREASE DESIRED                                

1 VARY O FIX   FITTED CONS  ERROR IN CONS                      

    0  0.5000000000000000       +- 0.000E+00

    1  -1.502885941371286       +- 0.652E-01

    1   1.851590712127221       +- 0.113    

    1  -1.350498684163228       +- 0.964E-01

    1  0.6470300251477153       +- 0.510E-01

    1 -0.2101536907945629       +- 0.178E-01

    1  0.4549469611716167E-01   +- 0.405E-02

    1 -0.6057434494797435E-02   +- 0.562E-03

    1  0.3854291140973252E-03   +- 0.371E-04

    1  -4.134151049838160       +- 0.130   

    1   7.368071342513580       +- 0.372   

    1  -7.633077301200458       +- 0.482   

    1   5.148920726637210       +- 0.372   

    1  -2.373642689043555       +- 0.188   

    1  0.7557095175197603       +- 0.640E-01

    1 -0.1619578118794217       +- 0.144E-01

    1  0.2147303973187547E-01   +- 0.199E-02

    1 -0.1366310760644399E-02   +- 0.131E-03

 CHIS= 4580.     FOR  5001 DATA POINTS

 CHI USED IN CALCULATING ERRORS IS  4984.   

 

The minimum without the truncation correction was 4519.  Simply re-reading the constants goes to 4619  4566 (This is the random/truncation zone). The constants above are the best for representing AiGauss.

Figure 4  Final fit to AiGauss.  -- Black is data, blue is fit, red is |data-fit|

Figure 5 (data-fit)/(data ×10-15)

            There is a truncation like error possibly in the Gauss Quadrature for -5.5<z<-3.5.  There is a small systematic effect at about -2.0 of a size < 1.  In general the “random” error in the fit is less than that in the data by a factor of (5000/18) = 16.7  that the error is < 10-15.

            The values below -6 are quickly generated to this accuracy using the asymptotic expansion ..\..\TestData\genint\Expanding the integrand.htm. doc. -- ..\..\TestData\genint\genint.zip Those above 0 are 1-the value below zero.  Thus the code to find AiGauss everywhere is ready to write.

Using the files

            The code for making this fit is in ..\..\zips\FittingBrack.zip

Directory of FittingBrack.zip

 aigau2.dir          -- sample direction file for fitting brack.out

AIGAUSS.OUT  -- data of the form exp(-z2)*Brack(z)  generated by direct integration of Gaussian

 BRACK.OUT  -- data of the form Brack(z)  Integral has been multiplied by exp(z2) with slight loss of accuracy

 

transmp.exe  -- Bailey’s code for converting to multi-precision -- ..\..\..\..\MultiplePrecision\bailey\Welcome.htm

mpfun.for  --- Bailey’s set of multi-precision routines

 

nlfit.wpj  Watcom ide’s for making nlfit

nlfit.tgt

cnlfit.bat  -- bat file for using wfl386 to make nlfit

 

mpnlfit.wpj  Watcom ide’s for making mpnlfit

mpnlfit.tgt

cmpnlfit.bat  bat file for using wfl386 to make mpnlfit

 

Fortran files

 MAOPEN.FOR  -- used by Watfor

MSOPEN.FOR  used by watcom

WSYSTEM.FOR  -- used by Watfor

openwsys.for  used by watcom

 

 nlfpade.for  - nlfit slightly modified to have Nnum, Nden, eps1, and eps2.  Passes flag for fitting Brack or exp*Brack to poly

MPnlfpade.for  Generated from nlfpade by running transmp.exe

 

 ROBMIN.FOR  --- minimization code.  Run transmp.exe and enter robmin.for to produce

 MProbmin.for  --- MProbmin

 

polyexpa.for  --- The poly routine for making exp(x2)×Brack(x) or simply Brack(x)