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.