Creating a Numerical Integration Function

[nova]

I'm a VERY basic user of Excel. I have a mathematics and programming background, but no knowledge or experience in VBA programming.

I'm wishing to create a new function as a VBA custom function procedure to numerically intergate an expression (such as 3*X^3 + 2*X^2), defined by a worksheet formula (not by data points). The five parameters required by the function would be in the form of cell references which point to (1) the expression to be integrated (such as the above example), (2) the integrating variable (in this example, X), (3) the lower limit of X, (4) the upper limit of X, and (5) a tolerance value for termination.

I've been trying to implement the Trapezoidal Method - but not a flicker of light so far!!.

I understand the Formula property of a cell may be used to get the expression to be integrated into the VBA code as text. Then the Substitute function may be used to replace the variable of interest by a value, and the Evaluate method may be used to get the new value of the expression. I have no experience in implementing these steps.

I've made extensive attempts to do the above, but without success so far. ("I started off with nothing, and still have most of it left"!).

I'm hoping someone may be able to assist me, or else point me to a relevant area of the forum (or elsewhere), by demonstrating the methods and techniques needed to achieve the function I'm seeking. Perhaps it has all been done before, but I haven't located it.

I'd be most grateful for any help.

 

[lalbatros]

nova,

Here is a quick and very dirty code that does roughly what you asked:

Function simpleInteg(fun, var, a, b)
    n = 100
    integ = 0
    For i = 0 To n
        v = a + (b - a) / n * i
        fun1 = Replace(fun, var, v)
        If i = 0 Then
            integ = integ + Evaluate(fun1) / 2
        ElseIf i = n Then
            integ = integ + Evaluate(fun1) / 2
        Else
            integ = integ + Evaluate(fun1)
        End If
    Next i
    simpleInteg = integ * (b - a) / n
End Function

Function cheatInteg(ifun, var, a, b)
    ifuna = Replace(ifun, var, a)
    ifunb = Replace(ifun, var, b)
    cheatInteg = Evaluate(ifunb) - Evaluate(ifuna)
End Function

Sub testExample()
    fun = "x"
    ifun = "x^2/2"
    a = 0
    b = 1
    Debug.Print simpleInteg(fun, "x", a, b), cheatInteg(ifun, "x", a, b)
    fun = "cos(x)"
    ifun = "sin(x)"
    a = 0
    b = 1
    Debug.Print simpleInteg(fun, "x", a, b), cheatInteg(ifun, "x", a, b)
    fun = "exp(x)"
    ifun = "exp(x)"
    a = 0
    b = 1
    Debug.Print simpleInteg(fun, "x", a, b), cheatInteg(ifun, "x", a, b)
End Sub

The last two functions are for testing purposes.

You will see that the last of the three examples will fail.
That is so because I simply subsitute the variables by using a string substitution.
This shows that to implement what you have requested, it is necessary to parse the expression in order to identify when the variable really appears in the expression. The evaluate function provided by excel has not ben designed to substitute variables and therefore substitution of variable is the difficult part of your question.

There is an excellent and free and open source solution to your question on this web page:

http://digilander.libero.it/foxes/mathparser/MathExpressionsParser.htm

I have been using this VBA expression parser quite intensively in some of my excel applications, and it is really excellent. The latest version is a bit slower but more powerful than previous versions.

 

[nova]

Thank you for your generous response and your effort to assist me. I now need to step through your solution and disect it line by line.

It appears my "simple" request is not so simple, ie quite an amount of code is needed. I have some work to do!!

Thank you again.

Kind regards.

 

[nova]

Hello lalbatros,<?xml:namespace prefix = o ns = "urn:schemas-microsoft-comfficeffice" /><o></o>
<o></o>
Thanks to your excellent advice, I have succeeded in producing a VB numerical integrator function for Excel. It also seems to be working ok!!<o></o>
<o></o>
The code, listed below, uses the Romberg method. However, I have added two refinements as follows: (1) it uses the Midpoint rule for function evaluations, and not the end points. This allows either or both limits of integration to be at singularities, etc. Very clever. (2) A non-linear transformation of X is included to achieve non-uniform spacing of sample points. This helps to prevent resonance phenomena occurring when periodic functions are integrated. Both of these refinements were implemented in the HP-34 calculator Integrate routine.<o></o>
<o></o>
I have added a fifth parameter for specifying a relative error value, such as .0000001. This is used for the convergence criteria test at the exit point in the code.<o></o>
<o></o>
For testing, I use the following (parameter) data:<o></o>
<o></o>
1. Expression = X^(1/8)<o></o>
2. Variable = X<o></o>
3. Lower limit a = 0<o></o>
4. Upper limit b = 1<o></o>
5. RelativeError = 0.0000001 (say)<o></o>
<o></o>
The true answer is 0.88888888888888888888.... etc, which allows for easy viewing and comparison of results.<o></o>
<o></o>
I'm hoping you may be able to advise me further regarding the program:<o></o>
<o></o>
Q1. Is there a way to improve the precision and accuracy of the results? As it stands, the code seems to be able to produce around 10 significant digits of accuracy. Is there a way to improve that to say 20 digits before round-off problems appear? I do not have knowledge of Excel's algorithms, accuracy aspects, or recommended calculating techniques.<o></o>
<o></o>
Q2. Are all the Dim and Type statements needed? Perhaps there are defaults for these. I note that your example did not use them, although other references seem to always show them being used.<o></o>
<o></o>
Q3. Are the Type statements correct for the respective uses of the variables, and for maximising speed, precision and accuracy?<o></o>
<o></o>
Q4. Are there, say, Excel environment aspects able to be considered and used to maximise the performance of the code?<o></o>
<o></o>
Q5. Are there other (coding) changes I could make to maximise the performance of the program?<o></o>
<o></o>
Q6. I'm wanting to be able to 'break-out' of the program if it seems to be in an endless loop, etc. Can the Esc key be programmed to do this? Are there other ways normally used to do this?<o></o>
<o></o>
I would be most grateful for any advice on the above points.<o></o>
<o></o>
Thank you again for your generous spirit. I have tried to put your advice to good use in the program, and it seems to be working very well indeed. However, this is also my first VB routine, so I expect it is at the kindergarten level, and much improvement awaits!!<o></o>
<o></o>
Kind regards,<o></o>
<o></o>
nova

<o:p></o:p>

<o:p></o:p>
<o:p></o:p>
===========================<o:p></o:p>
Start of Code<o:p></o:p>
===========================<o:p></o:p>
<o:p></o:p>
Function ROMBERG(Expression, Variable, a, b, RelativeError)<o:p></o:p>
Dim Sk As Double
Dim J As Integer, K As Integer
Dim TEMP1 As Double, DELTAU As Double, Ui As Double
Dim X As Double, Func As String, FX As Double
Dim M(30, 30) As Double<o:p></o:p>
<o:p></o:p>
Sk = 0
K = -1<o:p></o:p>
Do
Do
K = K + 1
TEMP1 = 2 ^ (-K)
DELTAU = 2 * TEMP1
Ui = -1 + TEMP1<o:p></o:p>
Do
X = (b - a) / 4 * Ui * (3 - Ui ^ 2) + (b + a) / 2<o:p></o:p>
Func = Replace(Expression, Variable, X)<o:p></o:p>
FX = Evaluate(Func)<o:p></o:p>
Sk = Sk + FX * (1 - Ui ^ 2)
Ui = Ui + DELTAU
Loop While Ui < 1
M(K, 0) = 3 * (b - a) / 4 * TEMP1 * Sk<o:p></o:p>
Loop While K = 0<o:p></o:p>
For J = 1 To K
M(K, J) = M(K, J - 1) + ((M(K, J - 1) - M(K - 1, J - 1)) / (4 ^ J - 1))
Next J<o:p></o:p>
Loop While Abs(M(K, K) - M(K - 1, K - 1)) / Abs(M(K, K)) > <o:p></o:p>
RelativeError<o:p></o:p>
ROMBERG = M(K, K)<o:p></o:p>
End Function<o:p></o:p>
<o:p></o:p>
=============================<o:p></o:p>
End of Code<o:p></o:p>
=============================<o:p></o:p>
<o:p>

</o>
<o></o>
<o>EDIT: tags added by mod</o>

 

[nova]

Sorry lalbatros, but all my code formatting just went out the window when I posted my reply. I need to look further into that so it doesn't happen again.

nova

 

[DonkeyOte]

nova, if you can edit your penultimate post...

type

[ code ]

before first line of VBA

and

[ /code ]

after last line... note: remove spaces in between [ ]

HTH

 

[nova]

Thank you for the advice. All understood.

nova

 

[lalbatros]

Hi Nova!

First of all I should recommend you "Numerical Recipes" (NR) which is one of the most used books on my desk.
See: http://www.nrbook.com/a/bookfpdf.php or: http://www.fizyka.umk.pl/nrbook/bookcpdf.html

Q1
You need first to find out the accuracy with the double type of visual basic.
You can do a simple test like calculating ((1+1E-15) - 1).
You will see that ((1+1E-16) - 1) returns zero, which gives you and idea of the precision.
All commercial numerical libraries need the values of some "machine constant" to run properly or optimally. See chapters 1.2 and 20.1 in NR.
The accuracy can be improved by several ways but there are always limitations.
For example, reducing the integration step will often improve the accuracy, but if the step is smaller than some limit it will lose any sense for many reasons (computation time, value smaller than the smallest number that can be represented ...)
However, it is possible to go very far in precision since it is only a matter of computer resource and programming.
One of the programming methods is called "Arbitray precision computing". You can read about that in NR chap 20.6. In this case, you don't rely anymore on pre-defined type of your programming language. A "multiprecision library" has been developped for excel, you can find it there: http://digilander.libero.it/foxes/SoftwareDownload.htm (same people as clsMathParser)

Q2
If you want a minimal level of wuality in your program, put this statement at the beginning of each module:

Option Explicit

And correct your code any time VB complains.
This will at least make thing clear in your program, but in VB this will also often dramatically accelerate your code. I have experienced that recently: by Typing each variables I improved significiantly the speed of my code (say x2), then using "VB Watch" I could find out many other optimisations and my code finally was running 10 times faster.

When I said dirty programming, this is FIRST OF ALL because nothing was dim-ed !!!
The second reason was that the code was not much re-usable, with little abstraction. this is typically what one does to get quick result, but such code is short-lived usually.

For you problem, I consider that writing an integration routine specically working for string-defined function is nonsense and dirty. The fast that your need is like that does not mean that you need to implement it literally in this way. And this is really the most dirty aspect.

Actually, you can find a lot of integration routine on the web. Some have been used since very long and improved. Writing such an algortihm makes sense only for an exercice or for creating a new algorithm. Therefore you should reuse an existing algorithm as much as possible. Integration routines are essentially numerical algorithm (but it could be symbolic too, like with Mathematica or Mapple or Maxima, ...). If you approach is numerical, then use tested numerical integration routines and design the interface you need.

This brings us to the topic of "how to pass a function as argument to a routine".
This was already solved by Fortran in the 60's.
Since then, the techniques have evolved enormously.
It is clear that today, a function must been seen as an object, define in VB in a Class module.
Typically a function is object that needs variable and that can at least return a value.
You can pass this object ot you integration routine.
The clsMathparser that I mentioned is precisely an implementation of a "function" based on string definitions.

You could already improve your code by "encapsulating" the "Evaluate" thing in your code.
This means removing this excel-specific stuff from you integration routine.
Put something more general in your integration routine.
Later, if you want to work with the "Evaluate" from excel creat a small class module for that.
But, if you like the clsMathParser I mentioned, you should be able to use it without changing your integration routine at all.
That would be an exercice about polymorphism, I think this is the name.

Q3
No sure about the meaning of your question.
Always type your variable.
If you do not, the compiler work will be delayed till run-time, and will even be much less efficient.

Q4
Of course!
You should never rely on any excel-things unless the goal is totally excel-thing.
An integration routine is totally unrelated to excel, therefore don't rely on anything from excel.

Of course, at some point, your stuff need to run within excel (this is requested), then you need to create the needed interface.
Again, the "Evaluate" function is an example: don't use it in your integration routine, instead use a standard way or a more general way of using/defining the function to integrate. Later, when using your coe within excel, add a small amount of code to work with the excel function "evaluate" if you need it or if you want it.

Sorry for being long and dirty.

Let me stress again, that I would personally never use the "Evaluate" function from Excel.
This is for the reasons I explained already, but also for other problem I encoutered with localisation. I did once an xl program that worked fine for my EU collegues but not anymore in the US! Simply because of the representaion of numbers. It made absolutely no sense to remedy that. Instead I went back to my real requirements. This is only one of numerous reasons why I finally opted for clsMathParser.

Another reason is the speed when the same function must be evaluated repeatedly, which is your case too.

Q5
Analyse the timing with "VB Watch" or with a self made tool (inserting time spies).
http://www.aivosto.com/vbwatch.html
Unfortunately, this is not available for Excel-VBA, I believe.
But if you own VB6, your could try most of your stuff in VB6. Excel is not really needed.

Q6
Many possibilities probably.
Maybe start a timer that will intervene in time.

 

[nova]

Hi lalbatros,

You've provided excellent advice and background information and detail. You are right - I do need to start on the general pathway with such routines, to keep it independent of Excel (or whatever). So much to learn - so little time!!

Thank you for your responses to my questions. In a short space of time, you have enabled me to program the basic functionality I was seeking. I now need to expand it and surround it with a universal framework, so it can be used "anywhere" with little extra effort.

Once again, you have been very generous with your time and effort to help me. It is nice to have a window into a group of people who care so much about newbies. Thank you.

I now need to locate a copy of the NR book. It seems to be a great resource.

Thank you once again.

Kind regards.

 

[lalbatros]

Nova,

I have put in a sample excel file what I explained in my last post.
You will find it there:

http://www.geocities.com/l.albatros/temp/integSample.xls

For a state-of-the-art approach, have a look at the Visual Numerics math package:

http://www.vni.com/products/imsl/jmsl/v50/manual/api/com/imsl/math/package-summary.html

You will see that functions are often represented as objects in this class library.
Often the numerical work to be performed not only needs numerical data but also function definitions. The "function object" contain the definition and can be accessed by the algorithms in a standardized way. The actual implementation of these object only need to satisfy the pre-defined interface requirement, but very different implementations can be defined and used as long as the predefined interface is available. For example such functions could be implemented totaly based on text strings (as in your question), but it could be based on an excel sheet (one of the examples I supplied) or it could be a formula in VB or it could be a complex engineering software properly interfaced.

Hope you enjoy the 21st century programming applied to numerical computations. ;>}