02-08-2017, 02:19 PM (This post was last modified: 02-08-2017 10:23 PM by Namir.)
Post: #12
 Namir Senior Member Posts: 1,012 Joined: Dec 2013
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.
 « Next Oldest | Next Newest »