The tanh-sinh quadrature - Printable Version

Code:



from mpmath import *

import time



# number of digits

mp.dps = 16

# repeat this # of times (to get better time estimations)

nloops = 100



# test equations

# equation =  0 => x*log(1 + x); 0; 1; 1/4

#             1 => x**2*atan(x); 0; 1; (pi - 2 + 2ln2)/12

#             2 => exp(x)cos(x); 0; pi/2; (exp(pi/2) - 1)/2

#             3 => atan(sqrt(2 + x**2))/(1 + x**2)/sqrt(2 + x**2); 0; 1; 5pi**2/96

#             4 => sqrt(x)ln(x); 0; 1; -4/9

#             5 => sqrt(1 - x**2); 0; 1; pi/4

#             6 => sqrt(x)/sqrt(1 - x**2); 0; 1; 2sqrt(pi)Gamma(3/4)/Gamma(1/4)

#             7 => ln(x)**2; 0; 1; 2

#             8 => ln(cos(x)); 0; pi/2; -pi*ln(2)/2

#             9 => sqrt(tan(x)); 0; pi/2; pi*sqrt(2)/2



print



# try all equations

for equation in range(10):

  if equation == 0:

    # limits

    a = mpf('0')

    b = mpf('1')

    # function

    def f(x):

      return x*ln(1 + x)

    # expected result (not used during computation)

    s = mpf('1/4')

  if equation == 1:

    a = mpf('0')

    b = mpf('1')

    def f(x):

      return x**2*atan(x)

    s = (pi() - 2 + 2*ln(mpf('2')))/12

  if equation == 2:

    a = mpf('0')

    b = pi()/2

    def f(x):

      return cos(x)*exp(x)

    s = (exp(pi()/2) - 1)/2

  if equation == 3:

    a = mpf('0')

    b = mpf('1')

    def f(x):

      return atan(sqrt(2 + x**2))/(1 + x**2)/sqrt(2 + x**2)

    s = 5*pi()**2/96

  if equation == 4:

    a = mpf('0')

    b = mpf('1')

    def f(x):

      return sqrt(x)*ln(x)

    s = mpf('-4/9')

  if equation == 5:

    a = mpf('0')

    b = mpf('1')

    def f(x):

      return sqrt(1 - x**2)

    s = pi()/4

  if equation == 6:

    a = mpf('0')

    b = mpf('1')

    def f(x):

      return sqrt(x)/sqrt(1 - x**2)

    s = 2*sqrt(pi())*gamma(mpf('3/4'))/gamma(mpf('1/4'))

  if equation == 7:

    a = mpf('0')

    b = mpf('1')

    def f(x):

      return ln(x)**2

    s = mpf('2')

  if equation == 8:

    a = mpf('0')

    b = pi()/2

    def f(x):

      return ln(cos(x))

    s = -pi()*ln(mpf('2'))/2

  if equation == 9:

    a = mpf('0')

    b = pi()/2

    def f(x):

      return sqrt(tan(x))

    s = pi()*sqrt(mpf('2'))/2



  # to measure algorithm execution time

  tt0 = time.time()

  # repeat nloops times to get better time estimations

  for m in range(nloops):

    # counters

    tnfe = 0  # counts function evaluations

    hyp = 0   # counts hyperbolics



    # we need a < b

    a, b = (a, b) if b > a else (b, a)

    # x = bpa2 + bma2*r

    bpa2 = (b + a)/2

    bma2 = (b - a)/2



    # epsilon

    eps = mpf('10')**(-mp.dps)

    # convergence threshold

    thr = mpf('10')**(-mp.dps/2)



    pi2 = pi()/2



    # (approx) maximum t that yields

    # the maximum representable r & x

    # values strictly below the upper

    # limits +1 (for r) and b (for x)

    tm = asinh(ln(2*min(1, bma2)/eps)/(2*pi2))

    hyp += 2



    # level

    k = 0

    # maximum allowed level

    maxlevel = int(ceil(log(mp.dps, 2))) + 2

    # ss is the final result

    # 1st addend at order 0

    ss = f(bpa2)

    tnfe += 1

    # "initial" previous computed result, used

    # in the convergence test

    sshp = 2*ss + 1

    # progress thru levels

    while k <= maxlevel:

      h = mpf('2')**-k

      N = int(floor(tm/h))

      j = 1

      while j <= N:

        t = j*h

        csh = pi2*sinh(t)

        ch = cosh(t)

        w = ch / cosh(csh)**2

        r = tanh(csh)

        hyp += 4

        fp = f(bpa2 + bma2*r)

        fm = f(bpa2 - bma2*r)

        p = w*(fp + fm)

        tnfe += 2

        ss += p

        # at level 0 must sweep all j values,

        # at other levels only the odd ones

        j += 2 if k > 0 else 1

      # converged?

      if abs(2*sshp/ss - 1) < thr: break

      # no, advance to next level

      k += 1

      # store the just computed approximation

      # for the next convergence test

      sshp = ss

    # apply constant coefficients

    ss *= bma2*pi2*h

    # done, measure time

    tt1 = time.time()



  # print results

  print equation, tnfe, hyp, k, (tt1 - tt0)/nloops, ss, ss - s, ss/s - 1

Code:

Option Explicit



' This program is based on the Python code shared by "Cesar" who goes by

' the user name emece67 on the hpmususem.com web site.



Function Fx(ByVal sFx As String, ByVal X As Double) As Double

  sFx = Replace(sFx, "EXP(", "!")

  sFx = Replace(sFx, "X", "(" & X & ")")

  sFx = Replace(sFx, "!", "EXP(")

  Fx = Evaluate(sFx)

End Function



Function Min(X, Y)

  Min = WorksheetFunction.Min(X, Y)

End Function



Function Floor(ByVal X As Double) As Double

  Floor = WorksheetFunction.Floor_Math(X)

End Function



Function Ceil(ByVal X As Double) As Double

  Ceil = WorksheetFunction.Ceiling_Math(X)

End Function



Function Sinh(ByVal X As Double) As Double

  Sinh = WorksheetFunction.Sinh(X)

End Function



Function Cosh(ByVal X As Double) As Double

  Cosh = WorksheetFunction.Cosh(X)

End Function



Function ASinh(ByVal X As Double) As Double

  ASinh = WorksheetFunction.ASinh(X)

End Function



Function Tanh(ByVal X As Double) As Double

  Tanh = WorksheetFunction.Tanh(X)

End Function



Sub WriteToColumnA()

' set up tags of columns A

  [A1].Value = "A"

  [A2].Value = "B"

  [A3].Value = "F(x)"

  [A4].Value = "Exact Result"

  [A5].Value = "Result"

  [A6].Value = "Err"

  [A7].Value = "%Err"

  [A8].Value = "Fx Calls"

  [A9].value = "Hyperbolics Count"

End Sub



Function TanhSinhQuadrature(ByVal sFx As String, ByVal A As Double, ByVal B As Double, _

              ByRef numFxCalls As Long, ByRef hyperbolicsCount As Long) As Double

  Dim K As Long, J As Long, Col As Integer

  Dim BuffAB As Double

  Dim ResultExact As Double, BmA2 As Double, BpA2 As Double, eps As Double

  Dim convergenceThreshold As Double, halfPi As Double, tMax As Double, h As Double

  Dim dps As Byte, maxLevel As Long, N As Long

  Dim ch As Double, csh As Double, w As Double, r As Double, fp As Double, fm As Double

  Dim p As Double, resultFinal As Double, resultPrevious As Double, t As Double

 

 

  If A > B Then

    BuffAB = A

    A = B

    B = BuffAB

  End If

  

  numFxCalls = 0

  hyperbolicsCount = 0

  dps = 16

  BpA2 = (B + A) / 2

  BmA2 = (B - A) / 2

  eps = 1 / 10 ^ dps

  ' convergence threshold

  convergenceThreshold = 1 / 10 ^ (dps / 2)

  halfPi = 2 * Atn(1)

  ' (approx) maximum t that yields

  ' the maximum representable r & x

  ' values strictly below the upper

  ' limits +1 (for r) and b (for x)

  tMax = ASinh(Log(2 * Min(1, BmA2) / eps) / (2 * halfPi))

  hyperbolicsCount = hyperbolicsCount + 2

  ' level

  K = 0

  ' maximum allowed level

  maxLevel = Int(Ceil(Log(dps) / Log(2))) + 2

  ' resultFinal is the final result

  ' 1st addend at order 0

  resultFinal = Fx(sFx, BpA2)

  numFxCalls = numFxCalls + 1

  ' "initial" previous computed result, used

  ' in the convergence test

  resultPrevious = 2 * resultFinal + 1

  ' progress thru levels

  Do While K <= maxLevel

    DoEvents

    h = 1 / 2 ^ K

    N = Int(Floor(tMax / h))

    J = 1

    Do While J <= N

      DoEvents

      t = J * h

      csh = halfPi * Sinh(t)

      ch = Cosh(t)

      w = ch / Cosh(csh) ^ 2

      r = Tanh(csh)

      hyperbolicsCount = hyperbolicsCount + 4

      fp = Fx(sFx, BpA2 + BmA2 * r)

      fm = Fx(sFx, BpA2 - BmA2 * r)

      p = w * (fp + fm)

      numFxCalls = numFxCalls + 2

      resultFinal = resultFinal + p

      ' at level 0 must sweep all j values,

      ' at other levels only the odd ones

      If K > 0 Then J = J + 2 Else J = J + 1 ' j += 2 if k > 0 else 1

    Loop

    ' converged?

    If Abs(2 * resultPrevious / resultFinal - 1) < convergenceThreshold Then Exit Do

    ' no, advance to next level

    K = K + 1

    ' store the just computed approximation

    ' for the next convergence test

    resultPrevious = resultFinal

  Loop

  ' apply constant coefficients

  resultFinal = BmA2 * halfPi * h * resultFinal

  TanhSinhQuadrature = resultFinal

End Function



Sub goColumns()

  Dim Col As Integer

  Dim A As Double, B As Double, sFx As String

  Dim ResultExact As Double

  Dim numFxCalls As Long, hyperbolicsCount As Long

  Dim resultFinal As Double

  

  ' you can comment the nest statement after you run the code once

  Call WriteToColumnA

  

  On Error GoTo HandleErr

  Col = 2

  Do Until IsEmpty(Cells(2, Col))

    A = Cells(1, Col)

    B = Cells(2, Col)

    sFx = UCase(Cells(3, Col))

    sFx = Replace(sFx, " ", "")

    ResultExact = Cells(4, Col)

    

    resultFinal = TanhSinhQuadrature(sFx, A, B, numFxCalls, hyperbolicsCount)



    Cells(5, Col) = resultFinal

    Cells(6, Col) = ResultExact - resultFinal

    Cells(7, Col) = 100 * Cells(6, Col) / ResultExact

    Cells(8, Col) = numFxCalls

    Cells(9, Col) = hyperbolicsCount

GoHere:

    Col = Col + 1

  Loop

    

ExitProc:

  Exit Sub

  

HandleErr:

  ' uncomment nest statement and set as breakpoint when you want to debug

  ' Resume Next

  Resume GoHere

    

End Sub