Option Explicit
' shg 2011
' Scope Name Description
' ------- ---------- -----------------------------------------------------
' Public Integral Returns the definite integral of a function over a
' specified interval; sets up the function used by Func.
' Uses ArgReplace & SimpsonInt
' Public FuncEval Evaluates a function with constant substitution. Uses
' ArgReplace.
' Private ArgReplace Replaces constants in an expression with values
' Private SimpsonInt Simpson integrator for a previously set up function
' Private Func Evaluates a previously setup function for a particular x
Dim oRE As Object
Dim sFunc As String
Function Func(x As Double) As Double
Debug.Print oRE.Replace(sFunc, x)
Func = Evaluate(oRE.Replace(sFunc, x))
End Function
Function Integral(Expression As String, _
xInitial As Double, xFinal As Double, Intervals As Long, _
ParamArray Constants() As Variant) As Variant
' shg 2011
' UDF or VBA
' Returns the numerical integration of Expression using the Composite
' Simpson's Rule
' Expression Expression to integrate. Any valid Excel real-number
' expression; must include the integration variable "x", and
' may optionally contain constants.
' xInitial Initial value of x
' xFinal Final value of x
' Intervals Number of intervals used for integration (will be rounded up
' to an even number). The value should be sufficiently large
' to achieve required accuracy. Larger values will take longer
' to compute.
' Constants Optional; if included, must be an alternating list of
' strings (constant names) and numbers (constant values), one
' pair for each constant in Expression.
' Examples =Integral("a*sin(x^2)+b*cos(x)", 0, 1, 100, "a", $A2, "b", B$1)
' =Integral("if(x=0, 1, sin(x)/x)", 0, PI(), 100)
Dim nConst As Long
Dim iConst As Long
Dim sNam As String
Dim dVal As Double
If Intervals < 1 Then GoTo BadnInt
If oRE Is Nothing Then
Set oRE = CreateObject("VBScript.Regexp")
oRE.Global = True
oRE.IgnoreCase = True
End If
With oRE
.Pattern = "\bx\b"
If Not .Test(Expression) Then GoTo NoX
If ArgReplace(Expression, CVar(Constants)) = False Then
Integral = Expression
Else
' ArgReplace changes the pattern, so set it again
.Pattern = "\bx\b"
sFunc = Expression
On Error GoTo BadFunc
Integral = SimpsonInt(xInitial, xFinal, (Intervals + 1) \ 2)
' See note in AdaptiveSI below
' Integral = AdaptiveSI(xInitial, xFinal)
End If
End With
Exit Function
BadnInt:
Integral = "Intervals must be > 0"
Exit Function
NoX:
Integral = "Integration variable ""x"" does not appear in Expression"
Exit Function
' Error handler for Func
BadFunc:
Integral = "Unable to evaluate " & sFunc
Exit Function
End Function
Public Function FuncEval(Expression As String, _
ParamArray Constants() As Variant) As Variant
' shg 2011
' UDF Only
ArgReplace Expression, CVar(Constants)
FuncEval = Evaluate(Expression)
If IsError(FuncEval) Then FuncEval = "Unable to evaluate " & sFunc
End Function
Private Function ArgReplace(Expression As String, _
Constants As Variant) As Boolean
' shg 2011
' VBA only
' Replaces the named arguments in Expression with values
' Expression Expression containing named constants
' Constants Alternating list of strings (constant names) and numbers
' (constant values), one pair for each constant in Expression.
' Returns True if parens in Expression are balanced and all substitutions made
' Otherwise returns False with error message in Expression.
Dim nConst As Long
Dim iConst As Long
Dim sNam As String
Dim dVal As Double
If Len(Replace(Expression, "(", "")) <> Len(Replace(Expression, ")", "")) Then GoTo BadParens
nConst = UBound(Constants) + 1
If nConst Mod 1 Then GoTo BadPramList
If oRE Is Nothing Then
Set oRE = CreateObject("VBScript.Regexp")
oRE.Global = True
oRE.IgnoreCase = True
End If
With oRE
For iConst = 0 To nConst - 1 Step 2
' parameters must be {string, number, string, number, ...}
If VarType(Constants(iConst)) <> vbString Or _
VarType(Constants(iConst + 1)) <> vbDouble Then
GoTo BadPramList
Else
sNam = Constants(iConst)
dVal = Constants(iConst + 1)
.Pattern = "\b" & sNam & "\b"
If Not .Test(Expression) Then GoTo NoConstant
Expression = .Replace(Expression, dVal)
End If
Next iConst
End With
ArgReplace = True
Exit Function
BadParens:
Expression = "Expression contains unbalanced parens"
Exit Function
BadPramList:
Expression = "Constants must be alternating strings (names) and numbers (values)"
Exit Function
NoConstant:
Expression = "Constant """ & sNam & """ does not appear in Expression"
Exit Function
End Function
'============= G e n e r a l S i m p s o n I n t e g r a t o r =============
Private Function SimpsonInt(xi As Double, xf As Double, ByVal n As Long) As Double
' shg 2006
' Returns the integral of Func (specific to problem) over the interval xi to xf
' using Composite Simpson's Rule over 2n intervals
Dim i As Double ' index
Dim dH As Double ' step size
Dim dF As Double ' cumulative integration of f
If n < 1 Then Exit Function
If xi = xf Then Exit Function
n = 2 * n
dH = (xf - xi) / n
Debug.Print xi, Func(xi)
dF = Func(xi) ' initial value
For i = 1 To n - 1 Step 2
Debug.Print xi + i * dH, Func(xi + i * dH)
dF = dF + 4 * Func(xi + i * dH) ' odd values
Next i
For i = 2 To n - 2 Step 2
Debug.Print xi + i * dH, Func(xi + i * dH)
dF = dF + 2 * Func(xi + i * dH) ' even values
Next i
Debug.Print xf, Func(xf)
dF = dF + Func(xf) ' final value
SimpsonInt = dF * dH / 3
End Function
Private Function AdaptiveSI(xi As Double, xf As Double) As Double
' shg 2011
' Returns the integral of Func (specific to problem) over the interval xi to xf
' using Composite Simpson's Rule.
' This works fine, other than trying to figure out how to limit the number of
' bisections. The interface to Integral would need to change to eliminate the
' Intervals argument.
Dim n As Long ' number of integration intervals
Dim dH As Double ' interval width
Dim i As Double ' index
' contributions of ...
Dim dInF As Double ' initial and final values of f
Dim dOdd As Double ' odd values of f
Dim dEvn As Double ' even values of f
Dim dDlt As Double ' delta in successive approximations
If xi = xf Then Exit Function
' calculate for two intervals
n = 2
dH = (xf - xi) / n
dInF = (Func(xi) + Func(xf)) * dH / 3
dOdd = 4 * Func(xi + dH) * dH / 3
AdaptiveSI = dInF + dOdd
' Debug.Print AdaptiveSI
' bisect repeatedly
Do
AdaptiveSI = AdaptiveSI + dDlt
n = n * 2
dH = (xf - xi) / n
dInF = dInF / 2
dEvn = dOdd / 4 + dEvn / 2
dOdd = 0
For i = 1 To n - 1 Step 2
dOdd = dOdd + 4 * Func(xi + i * dH) * dH / 3 ' odd values
Next i
dDlt = dInF + dOdd + dEvn - AdaptiveSI
' Debug.Print n, dDlt, AdaptiveSI
Loop Until AdaptiveSI = AdaptiveSI + dDlt Or n >= 1024
End Function
