# 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,