Numerical Methods and Excel VBA

Posted in hacks and kludges on Thursday, October 08 2015

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 Ctrl-Shift-Enter).

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 fanno to 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.

tags: excel, vba,



comments powered by Disqus