The tanh-sinh quadrature

[Namir]

The following Excel VBA code is based on emece67's Python code.

The VBA code writes the tags for column A. The code uses columns B and up to solve one or more numerical integration problems. The input rows in column B and up contain the following data:

1. Row 1 contains the value of A, the lower integration limit.
2. Row 2 contains the value of B, the upper integration limit.
3. Row 3 contains text that describes the function to be integrated. For example the input can be "x*ln(1 + x)
". The input is case insensitive.
4. Row 4 contains the value of the exact integral. You can enter either a constant or a calculated value, such as "=LN(10)" or "=LN(B2/B1)".

Note: Make sure that math functions appearing in row 3 are properly interpreted by VBA's runtime system. You make need to change the function's name from, say, LN to LOG to avoid runtime error.

The output rows contain the following data:

1. Row 5 contains the calculated integral.
2. Row 6 contains the error in the integral.
3. Row 7 contains the percent error in the integral.
4. Row 8 contains the number of function calls.
5. Row 9 contains the number of "hyperbolic count" which is (usually?) twice the number of function call.

Here is the VBA code.

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

The above code should be a good step into generating HP Prime code for the tanh-sinh quadrature.