I'm a big fan of using circular references and other tricks to get excel to solve equations numerically and automatically. Nobody likes constantly clicking through the solver or using goal seek when the spreadsheet could just solve itself. However circular references can be a hassle, especially for the uninitiated.
Recently I've been taking some exel files I made for doing compressible flow calculations and mushing them together into a template with macros to do all the scut work. Compressible flow in pipes typically requires numerically solving equations and to automate this I used circular references and Newton's method. This works OK if the initial guess is pretty good, and if the user knows how to re-start calculations when things go awry (and all the cells turn into
#Value errors or something). But I want something that is more robust and more idiot proof.
What I would like to do is use Brent's method but as a user defined function, so someone can type
=root_of_equation(blah, blah, blah) and just get the answer. Turns out this isn't as easy as it sounds!
VBA has limitations, especially if you are used to languages with sensible class inheritance like Python, or where you can pass functions by reference like Python, or if you just like Python (I like Python). In Python I would so something like this:
def brent(f, lowerbound, upperbound): ''' implement brent's method here ''' def fanno(Ma, *args): ''' implement the compressible flow equation that needs to be solved here ''' solution = brent(lambda x: fanno(x, parameters), lowerbound, upperbound)
In VBA I have do something a little more round about. Firstly I create a class module
Brent for Brent's method with the method
solve, which takes an object
f and a set of bounds and finds the root of
f.eval(x). This is a round-about way of effectively passing one function to another function, except in VBA land I can only pass objects to other objects. One thing to note is that
brent.solve() returns a size 2 array containing the solution and the relative error of the solution.
Option Explicit Private pMAX_ITER As Integer Private pERROR_TOL As Double Public Property Get MAX_ITER() As Integer MAX_ITER = pMAX_ITER End Property Public Property Let MAX_ITER(value As Integer) pMAX_ITER = value End Property Public Property Get ERROR_TOL() As Double ERROR_TOL = pERROR_TOL End Property Public Property Let ERROR_TOL(value As Double) pERROR_TOL = value End Property Public Function solve(f As Variant, lowguess As Double, highguess As Double) As Variant Dim a As Double Dim b As Double Dim c As Double Dim d As Double Dim fa As Double Dim fb As Double Dim fc As Double a = lowguess b = highguess fa = f.eval(a) fb = f.eval(b) 'Throw a Value error if the root is not within the bounds If fa * fb >= 0 Then solve = CVErr(xlErrValue) Exit Function End If 'If the roots are backwards, swap them If Math.Abs(fa) < Math.Abs(fb) Then Dim swap As Double swap = a a = b b = swap swap = fa fa = fb fb = swap End If c = a fc = f.eval(c) Dim mflag As Boolean mflag = True Dim i As Integer Dim err As Double i = 0 err = 1# Do While (i < MAX_ITER) And (err > ERROR_TOL) Dim s As Double Dim fs As Double If fa <> fc And fb <> fc Then s = ((a * fb * fc) / ((fa - fb) * (fa - fc))) + ((b * fa * fc) / ((fb - fa) * (fb - fc))) + ((c * fa * fb) / ((fc - fa) * (fc - fb))) Else s = b - fb * ((b - a) / (fb - fa)) End If 'Use bisection to find s If ((s >= (3 * a + b) / 4 And s <= b) Or (s <= (3 * a + b) / 4 And s >= b)) Or _ (mflag And Math.Abs(s - b) >= Math.Abs(b - c) / 2) Or _ (Not mflag And Math.Abs(s - b) >= Math.Abs(c - d) / 2) Or _ (mflag And Math.Abs(b - c) < ERROR_TOL) Or _ (Not mflag And Math.Abs(c - d) < ERROR_TOL) Then s = (a + b) / 2 mflag = True Else mflag = False End If fs = f.eval(s) d = c c = b If fa * fs < 0 Then b = s fb = f.eval(b) Else a = s fa = f.eval(a) End If If Math.Abs(fa) < Math.Abs(fb) Then swap = a a = b b = swap swap = fa fa = fb fb = swap End If i = i + 1 err = Math.Abs((b - a) / a) If fb = 0# Then err = 0 Loop Dim arr(1 To 2) As Double arr(1) = b arr(2) = err solve = arr End Function
I also created a class module
Fanno for, in this example, solving fanno flow. The module implements the
eval function, so an instance of the
Brent class can take it and find its roots.
Option Explicit Private pk As Double Private pfp As Double Public Property Get k() As Double k = pk End Property Public Property Let k(value As Double) pk = value End Property Public Property Get fp() As Double fp = pfp End Property Public Property Let fp(value As Double) pfp = value End Property Public Function parameter(MachNumber As Double) As Double Dim Ma As Double Ma = MachNumber parameter = ((1 - Ma ^ 2) / (k * Ma ^ 2)) + ((k + 1) / (2 * k)) * Math.Log(((k + 1) * Ma ^ 2) / (2 + (k + 1) * Ma ^ 2)) End Function Public Function eval(MachNumber As Double) As Double eval = fp - parameter(MachNumber) End Function
At this point I have classes defining objects, but no code connecting the two and making it accessible as a user defined function in excel. To that end I wrote a method,
fanno_roots that lets me numerically solve the fanno flow equation for the Mach number using Brent's method.
Public Function fanno_root(fanno_parameter As Double, k As Double, lowbound As Double, highbound As Double) As Variant 'Initialize the brent method solver Dim brentq As Object Set brentq = New Brent brentq.MAX_ITER = 100 brentq.ERROR_TOL = 0.000000000001 'Initialize the fanno flow object Dim f As Object Set f = New Fanno f.k = k f.fp = fanno_parameter 'Calculate the result Dim result As Variant result = brentq.solve(f, lowbound, highbound) 'Set the resulting array to be the right direction (row or column) If Application.Caller.Rows.Count > 1 Then fanno_root = Application.Transpose(result) Else fanno_root = result End If End Function
This is an array function it will give you back the root of the equation and the relative error of the root if you give it the space (select 2 cells, enter
=fanno_root(fanno_parameter, k, lowguess, highguess) and then hit
I basically did the same thing in VBA as I did in Python, except that in Python functions are first class objects so I did not need to define a class for
brent and a class for
fanno, initialize them as objects, then pass
brent, I can just do it directly. The reason I went through all this hoop jumping is that I want my implementation of Brent's method to be portable and reusable, I have several points in this template (and other templates) where I want to use it.
The other big part of making this template more user friendly (rewinding back to my motivation) is automating the estimating of the boundaries for the root finding. This can be done with some understanding of the physics of the problem, but that is really a question for another time.