  Rank: Advanced Member Groups: Registered
Joined: 01/08/2010(UTC) Posts: 115   Location: Rome Was thanked: 13 time(s) in 13 post(s)
|
Hello Andrey In this period I have collected and studied the problem of numerical integration and I have seen that it is really a wide area of math in which is easy to get lost. So I have for now limited my purpose to the numerical integration of functions of one variable using non adaptative algorithms. I have found some solutions already implemented using Gauss quadrature method and have tested it in .....Excel that is the development environment I am more familiar with. Code:
Function GaussInt(Func As String, ParamA As Variant, ValueA As Variant, IntA As Variant)
Dim TOL As Double, Eform As String, LLim As Double, ULim As Double
Dim Result(1 To 1, 1 To 4) As Double, IERR As Long, ErrVal As Double
Dim XA() As Double, WA() As Double, EStep As String, Nnodes As Variant
Dim OldArea As Double, Area As Double, MaxLoops As Long
Dim I As Long, J As Long, K As Long, N As Long, IntAUB As Long
Dim X_1 As Double, x_3 As Double, X_2 As Double, W_2 As Double, SArea As Double
Dim W As Double
Result(1, 4) = MicroTimer
GetArray ParamA
GetArray ValueA
GetArray IntA
IntAUB = UBound(IntA) - LBound(IntA) + 1
If IntAUB < 3 Then
GaussInt = "Invalid IntA"
Exit Function
End If
Eform = Func
For I = 1 To UBound(ParamA, 1) - LBound(ParamA, 1) + 1
Eform = Replace(Eform, ParamA(I, 1), ValueA(I, 1))
Next I
LLim = IntA(2, 1)
ULim = IntA(3, 1)
If IntAUB > 3 Then TOL = IntA(4, 1) Else TOL = 0.0000000001
If IntAUB > 4 Then MaxLoops = IntA(5, 1)
If IsEmpty(MaxLoops) = True Or MaxLoops < 2 Then MaxLoops = 10
If IntAUB > 5 Then Nnodes = IntA(6, 1) Else Nnodes = 12
If Nnodes = "" Or Nnodes < 2 Or Nnodes > 12 Then Nnodes = 12
ReDim XA(1 To Nnodes)
ReDim WA(1 To Nnodes)
IERR = GetGaussA(Nnodes, XA(), WA())
OldArea = 0
N = 1
For K = 1 To MaxLoops ' Double the number of steps up to MaxLoops times
Area = 0
W = (ULim - LLim) / N
For I = 1 To N
X_1 = LLim + (I - 1) * W
x_3 = X_1 + W
X_2 = (x_3 + X_1) / 2
W_2 = (x_3 - X_1) / 2
For J = 1 To Nnodes
EStep = Replace(Eform, IntA(1, 1), X_2 + W_2 * XA(J))
SArea = Evaluate(EStep)
Area = Area + WA(J) * SArea
Next J
Next I
Area = Area * W_2
ErrVal = Abs((Area - OldArea) / Area)
If ErrVal < TOL Then GoTo Finish
OldArea = Area
N = 2 * N
Next K
K = K - 1
Finish:
Result(1, 1) = Area
Result(1, 2) = ErrVal
Result(1, 3) = K
Result(1, 4) = MicroTimer - Result(1, 4)
GaussInt = Result
End Function
Public Sub GetArray(arrayname)
Dim TEMP As Variant
If TypeName(arrayname) = "Range" Then
If arrayname.Rows.Count = 1 Then
If arrayname.Columns.Count = 1 Then
arrayname = Array(arrayname.Value2)
TEMP = arrayname(0)
ReDim arrayname(1 To 1, 1 To 1)
arrayname(1, 1) = TEMP
Else
arrayname = Array(arrayname.Value2)
arrayname = WorksheetFunction.Transpose(arrayname)
End If
Else
arrayname = arrayname.Value2
End If
ElseIf Not IsArray(arrayname) Then
arrayname = Array(arrayname)
arrayname = WorksheetFunction.Transpose(arrayname)
Else
arrayname = WorksheetFunction.Transpose(arrayname)
End If
End Sub
Function GetGaussA(Nnodes, ByRef XA, ByRef WA) As Long
Select Case Nnodes
Case 2
XA(1) = -0.577350269189626
XA(2) = 0.577350269189626
WA(1) = 1
WA(2) = 1
Case 3
XA(1) = -0.774596669241483
XA(2) = 0
XA(3) = 0.774596669241483
WA(1) = 0.555555555555555
WA(2) = 0.888888888888889
WA(3) = 0.555555555555555
Case 4
XA(1) = -0.861136311594053
XA(2) = -0.339981043584856
XA(3) = 0.339981043584856
XA(4) = 0.861136311594053
WA(1) = 0.347854845137453
WA(2) = 0.652145154862546
WA(3) = 0.652145154862546
WA(4) = 0.347854845137453
Case 5
XA(1) = -0.906179845938664
XA(2) = -0.538469310105683
XA(3) = 0
XA(4) = 0.538469310105683
XA(5) = 0.906179845938664
WA(1) = 0.236926885056189
WA(2) = 0.478628670499366
WA(3) = 0.568888888888889
WA(4) = 0.478628670499366
WA(5) = 0.236926885056189
Case 6
XA(1) = -0.932469514203152
XA(2) = -0.661209386466265
XA(3) = -0.238619186083197
XA(4) = 0.238619186083197
XA(5) = 0.661209386466265
XA(6) = 0.932469514203152
WA(1) = 0.171324492379171
WA(2) = 0.360761573048138
WA(3) = 0.467913934572691
WA(4) = 0.467913934572691
WA(5) = 0.360761573048138
WA(6) = 0.171324492379171
Case 7
XA(1) = -0.949107912342759
XA(2) = -0.741531185599395
XA(3) = -0.405845151377397
XA(4) = 0
XA(5) = 0.405845151377397
XA(6) = 0.741531185599395
XA(7) = 0.949107912342759
WA(1) = 0.129484966168869
WA(2) = 0.279705391489278
WA(3) = 0.381830050505117
WA(4) = 0.417959183673472
WA(5) = 0.381830050505117
WA(6) = 0.279705391489278
WA(7) = 0.129484966168869
Case 8
XA(1) = -0.960289856497536
XA(2) = -0.796666477413627
XA(3) = -0.525532409916329
XA(4) = -0.18343464249565
XA(5) = 0.18343464249565
XA(6) = 0.525532409916329
XA(7) = 0.796666477413627
XA(8) = 0.960289856497536
WA(1) = 0.101228536290376
WA(2) = 0.222381034453374
WA(3) = 0.313706645877887
WA(4) = 0.362683783378361
WA(5) = 0.362683783378361
WA(6) = 0.313706645877887
WA(7) = 0.222381034453374
WA(8) = 0.101228536290376
Case 9
XA(1) = -0.968160239507626
XA(2) = -0.836031107326637
XA(3) = -0.61337143270059
XA(4) = -0.324253423403808
XA(5) = 0
XA(6) = 0.324253423403808
XA(7) = 0.61337143270059
XA(8) = 0.836031107326635
XA(9) = 0.968160239507627
WA(1) = 8.12743883615759E-02
WA(2) = 0.180648160694854
WA(3) = 0.260610696402935
WA(4) = 0.312347077040002
WA(5) = 0.330239355001259
WA(6) = 0.312347077040002
WA(7) = 0.260610696402935
WA(8) = 0.180648160694857
WA(9) = 8.12743883615721E-02
Case 10
XA(1) = -0.973906528517172
XA(2) = -0.865063366688984
XA(3) = -0.679409568299024
XA(4) = -0.433395394129247
XA(5) = -0.148874338981631
XA(6) = 0.148874338981631
XA(7) = 0.433395394129247
XA(8) = 0.679409568299024
XA(9) = 0.865063366688984
XA(10) = 0.973906528517172
WA(1) = 6.66713443086868E-02
WA(2) = 0.149451349150573
WA(3) = 0.219086362515983
WA(4) = 0.269266719309996
WA(5) = 0.295524224714752
WA(6) = 0.295524224714752
WA(7) = 0.269266719309996
WA(8) = 0.219086362515983
WA(9) = 0.149451349150573
WA(10) = 6.66713443086868E-02
Case 11
XA(1) = -0.97822865814604
XA(2) = -0.88706259976812
XA(3) = -0.730152005574042
XA(4) = -0.519096129206811
XA(5) = -0.269543155952344
XA(6) = 0
XA(7) = 0.269543155952344
XA(8) = 0.519096129206811
XA(9) = 0.73015200557405
XA(10) = 0.887062599768093
XA(11) = 0.978228658146058
WA(1) = 5.56685671162158E-02
WA(2) = 0.125580369464874
WA(3) = 0.18629021092774
WA(4) = 0.233193764591992
WA(5) = 0.262804544510246
WA(6) = 0.2729250867779
WA(7) = 0.262804544510246
WA(8) = 0.233193764591993
WA(9) = 0.186290210927733
WA(10) = 0.125580369464913
WA(11) = 5.56685671161695E-02
Case 12
XA(1) = -0.981560634246732
XA(2) = -0.904117256370452
XA(3) = -0.769902674194317
XA(4) = -0.587317954286614
XA(5) = -0.36783149899818
XA(6) = -0.125233408511468
XA(7) = 0.125233408511468
XA(8) = 0.36783149899818
XA(9) = 0.587317954286614
XA(10) = 0.769902674194317
XA(11) = 0.904117256370452
XA(12) = 0.981560634246732
WA(1) = 4.71753363864754E-02
WA(2) = 0.106939325995363
WA(3) = 0.160078328543358
WA(4) = 0.203167426723067
WA(5) = 0.233492536538353
WA(6) = 0.249147045813402
WA(7) = 0.249147045813402
WA(8) = 0.233492536538353
WA(9) = 0.203167426723067
WA(10) = 0.160078328543358
WA(11) = 0.106939325995363
WA(12) = 4.71753363864754E-02
End Select
End Function
The problem now it trying to traslate this functions in VB Express programs ; class module as explained in your tutorial. I have read it and downloaded the VB Express IDE but I am a little bit intimidated by this new development environment. So I have a little bit of questions to ask: In developing a numerical integration program is necessary to input same informations: 1-The algebraic expression of the function to integrate ( a text string) 2-The variable of integration ( a text string ) 3-The integration extremes ( numbers ) now in Excel I have cells where to load these informations as text string and then internally using the trick of the evaluate method I can for different values of the variable of integration calculate the corresponding values of the function to integrate. What in my ignorance I am not able to reproduce ( for now ) in VB express is this first steps. In particular the one of analyzing the string that describes the Integral function and turn it in something that can be computed changing the values of the variable of integration. What I wold like to know is if I need to develop a math parser for this task or is something - as I guess - that is already implemented in Smath and so I need only to understand how I can access this internal description and how to use it for the number crunching algorithm. Because you have already developed a definite integral class/routine would be nice - if possible - to see how this conversion is implemented or how I can access the results of this conversion to use for further calculations. If is not possible to have these info's/part of code it would be enough to have just the header of the actual routine and the returning part so that I can focuse my attention to the core part that is related to the numerical integration only. Moreover I am not a class oriented but a more brutal stright to the result VBA programmer and this is creating me a quite steep learning curve in this new development environment. Sorry again for my newbie questions and best regards Francesco
|
