bromberg method compared with tanh-sinh and gauss-legendre methods



On Nov. 29, 2012, I wrote:
------------------------------------------------
>Comparison of bromberg method with the other
>available bigfloat methods (see Ch.9)
>for bigfloat quadrature of finite
>domain integrals.
----------------------------------------
The previous report concentrated on times required for
two test functions and for each function three values
of the working precision wp = local value of fpprec.

Some indications of the accuracy were implicitly
available in the printouts of the (rp = 20, wp = 40) runs.

Here is a summary of the accuracy for the two functions
chosen. Again, ts = tanh-sinh method,
gs = gauss-legendre method, and br = bromberg method.
-------------------------------------------------------------
  1.) well behaved monotonic function
            1/(1 + sqrt(x) over [1/2, 1]
      
 ===================================================
    accuracy comparisons (diff from true value)
 
  wp             ts                  gs                     br
     
  40       10^(-41)             0                   10^(-25)
  
  60       10^(-60)         10^(-60)          10^(-47)
  
  110      10^(-110)      10^(-109)        10^(-96)
  
=======================================================
  2.)  pure sinusoid function sin(50*x) over [0, 1]
  
===================================================
    accuracy comparisons (diff from true value)
 
  wp            ts                    gs                     br
 
  40       10^(-41)          10^(-40)         10^(-28)
  
  60       10^(-62)          10^(-59)         10^(-48)
  
  110      10^(-112)      10^(-109)        10^(-93)
  
=======================================================

In general, the bigfloat tanh-sinh and gauss-legendre
methods are more accurate than bromberg for a given wp.

------------------------------------------
For completeness, I show the coding for
mcompare.mac which recovers the times and
accuracy.
====================================
/* mcompare.mac

   functions:   
   br_test, ts_test, gs_test,
   mcompare2, mcompare3

   Each method-test prints a line
       method-code  abs(mval - tval)  #sec,
   and uses global tval = approx 'true value'
   to compute the absolute error.
   
 */  
   
load("c:/work2/nint.mac")$
load("c:/work2/tsquad.mac")$
load("c:/work2/gsquad.mac")$
load("c:/work2/qbromberg.mac")$

timer(tsquad_r, qbromberg, gsquad)$

/* bromberg method:
       br_test(e,x,x1,x2,rp,wp)
        calls qbromberg, which expects
        a real integrand.
 */

br_test (bee, bxx, bxx1, bxx2, brrp, bwwp):=
block ([mval, oldt, mt, mdiff],
    oldt : get (qbromberg,  'runtime),
    mval : qbromberg (bee, bxx, bxx1, bxx2, brrp, bwwp),
    mt : (get (qbromberg, 'runtime) - oldt)/100.0,
    if mval = false then mdiff:false
    else mdiff : block ([fpprec:bwwp], abs(mval - tval)),
    print(" br:  ", float(mdiff), "   ", mt))$
    
/* tanh-sinh method
           ts_test(e,x,x1,x2,rp,wp)
           calls tsquad_r which expects a
           real integrand.
 */
    
ts_test (tsee, tsxx, tsxx1, tsxx2, tsrrp, tswwp):=
block ([mval, oldt, mt, mdiff],
    oldt : get (tsquad_r, 'runtime),
    mval : tsquad_r (tsee, tsxx, tsxx1, tsxx2, tsrrp, tswwp),
    mt : (get (tsquad_r, 'runtime) - oldt)/100.0,
    if mval = false then mdiff:false
    else mdiff : block ([fpprec:tswwp], abs(mval - tval)),
    print(" ts:  ", float (mdiff), "  ", mt))$
    
    
/* gauss-legendre method
         gs_test(e,x,x1,x2,rp,wp)
         calls gsquad, which expects
         a real integrand.
 */
    
gs_test (gsee, gsxx, gsxx1, gsxx2, gsrrp, gswwp):=
block ([mval, oldt, mt, mdiff],
    oldt : get (gsquad, 'runtime),
    mval : gsquad (gsee, gsxx, gsxx1, gsxx2, gsrrp, gswwp),
    mt : (get (gsquad, 'runtime) - oldt)/100.0,
    if mval = false then mdiff:false
    else mdiff : block ([fpprec:gswwp], abs(mval - tval)),
    print(" gs:  ", float (mdiff), "  ", mt))$
    
    
/*     mcompare2 (e,x,x1,x2,rp,wp)
           calls ts_test and gs_test
 */
           
mcompare2 (ce,cx,cx1,cx2,crp,cwp) :=
        (ts_test (ce, cx, cx1, cx2, crp, cwp),
               gs_test (ce, cx, cx1, cx2, crp, cwp))$
      
/*     mcompare3 (e,x,x1,x2,rp,wp)
           calls mcompare2 and br_test.
 */
      
mcompare3 (cce, ccx, ccx1, ccx2, ccrp, ccwp) :=
     ( mcompare2 (cce, ccx, ccx1, ccx2, ccrp, ccwp),
       br_test (cce, ccx, ccx1, ccx2, ccrp, ccwp))$
==================================================
Ted Woollett