Non linear equations solver using VBA 1

[mwalshe]

I need to solve a number of non linear equations using VBA (specifically this is for estimation of the equilibrium composition of an gas mixture using Gibbs free energy minimization which will have at least 10 non lin. equations to solve). The problem is that the code I have (from "EXCEL FOR SCIENTISTS & ENGINEERS by E.Joseph Billo) using gaussian jordan elimination/matrix pivoting, it produces a zero pivot & thus an error when the matrix elements are divided by the pivot.

Can anyone suggest how I augment the code to overcome this problem (division by zero pivot term) or has anyone else a better code for non linear equations solver. CODE BELOW

''''''''''''''''''''CODE'''''''''''''''''''''''''''''''
Option Explicit
Option Base 1

Const ConvergenceTolerance = 0.00000001
Const IncermentNumericalDifferentiation = 0.000000001
Const Iterations = 50
Const expon = 2.718281828



Function SimultEqNL(equations, Variables, constants)
' Newton iteration method to find roots of nonlinear simultaneous equations
' Example:
' w^3 + 2w^2 + 3w + 4 = +12.828
' w.x + x.y + y.z = -3.919
' w^2 + 2w.x + x^2 = +1
' w + x + y - z = -3.663
'
' WHERE: constants = [12.828,-3.919,1,-3.663];

'On Error Resume Next

Dim i As Integer, j As Integer, k As Integer, N As Integer
Dim Nlterations As Integer
Dim R As Integer, C As Integer
Dim VarAddr() As String, FormulaString() As String
Dim con() As Double, A() As Double, B() As Double
Dim V() As Double
Dim Y1 As Double, Y2 As Double
Dim tolerance As Double, incr As Double

N = equations.Rows.count
k = Variables.Rows.count

If k = 1 Then k = Variables.Columns.count
If k <> N Then SimultEqNL = CVErr(xlErrRef): Exit Function
' Use the CVErr function to create user-defined errors in user-created procedures.
' For example, if you create a function that accepts several arguments and normally returns a string,
' you can have your function evaluate the input arguments to ensure they are within acceptable range.
' If they are not, it is likely your function will not return what you expect.
' In this event, CVErr allows you to return an error number that tells you what action to take.


ReDim VarAddr(N), FormulaString(N), V(N), con(N)
ReDim A(N, N + 1), B(N, N + 1)

tolerance = ConvergenceTolerance 'Convergence criterion.
incr = IncermentNumericalDifferentiation 'Increment for numerical differentiation.
Nlterations = Iterations

For i = 1 To N
VarAddr(i) = Variables(i).Address
' i.e. VarAddr(1) = $A$11
Next i

' Initial values
For i = 1 To N
con(i) = constants(i).Value
' Put the initial guesses into vector V()
V(i) = Variables(i).Value: If V(i) = 0 Then V(i) = 1
Next i

For j = 1 To Nlterations
' Create N x N matrix of partial derivatives.
For R = 1 To N ' n = equations.Rows.count
For C = 1 To N
' Formulastring is formula in which all but one variable in each equation is replaced by current values.
FormulaString(R) = Application.ConvertFormula(equations(R).Formula, xlA1, xlA1, xlAbsolute)
' xlA1 = Use xlA1 to return an A1-style reference. xlR1C1 = Use xlR1C1 to return an R1C1-style reference
' xlAbsolute = Convert to absolute row and column style

' ConvertFormula method used to convert cell references from A1 reference style to R1C1 reference style
For i = 1 To N
' Debug.Print FormulaString(R)
If i <> C Then FormulaString(R) = Application.Substitute(FormulaString(R), VarAddr(i), V(i))
' Substitutes new_text for old_text in a text string.
' FormulaString(R) = the reference to a cell containing text for which you want to substitute characters.
' Replace the address reference with the actual variable value i.e. $B$5^3 + 2*$B$5^2 + 3*$B$5 + 4-$R$5+3-1^3
' i = C means its on the diagonal of the matrix
Next i
' V() = vector of current variable values
If IsError(Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 + incr)))) Then MsgBox "ERROR IS FORMAULA EVALUATION"
If IsError(Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 - incr)))) Then MsgBox "ERROR IS FORMAULA EVALUATION"

' Y2 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), (V(C) + incr)))
' ' value of the equation at the current variable figure i.e. instead of evaluating the equation at VarAddr(C) evaluate at V(C)*(1+incr)
' Y1 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), (V(C) - incr)))
' A(R, C) = (Y2 - Y1) / (2 * incr)
Y2 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 + incr)))
' value of the equation at the current variable figure i.e. instead of evaluating the equation at VarAddr(C) evaluate at V(C)*(1+incr)
Y1 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 - incr)))
A(R, C) = (Y2 - Y1) / (2 * incr * V(C))
'
' Derivatives i.e Taylor Series approx (numerical differentiation) central difference:
' y' = (y,i+1 - y,i-1) / 2h where y' = derivative of y, h = step size
' F(x+h) = F(x) + hF'(x) thus F'(x) = [ F(x+h)-F(x) ] / h
Next C
Next R

'Augment matrix of derivatives with vector of constants.
For R = 1 To N
FormulaString(R) = Application.ConvertFormula(equations(R).Formula, xlA1, xlA1, xlAbsolute)
For C = 1 To N
FormulaString(R) = Application.Substitute(FormulaString(R), VarAddr(C), V(C))
Next C
If IsError(Evaluate(FormulaString(R))) Then MsgBox "ERROR IN FORMULA FO Augment matrix of derivatives"
A(R, N + 1) = con(R) - Evaluate(FormulaString(R))
Next R

For i = 1 To N
If Abs((A(i, N + 1)) / V(i)) > tolerance Then GoTo Refine
Next i

SimultEqNL = Application.Transpose(V)
Exit Function

Refine: Call GaussJordan3(N, A, B)
'Update V values
For i = 1 To N
V(i) = V(i) + A(i, N + 1)
Next i
'Debug.Print j, "", V(1), V(2), V(3), V(4) ', V(5), V(6), V(7), V(8), V(9)
Next j
' Exit here if no convergence after a specified number of iteration
SimultEqNL = CVErr(xlErrNA)

End Function
'. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sub GaussJordan3(N, AugMatrix, TempMatrix)
Dim i As Integer, j As Integer, k As Integer, L As Integer, LL As Integer, P As Integer, M As Integer, MM As Integer, MMM As Integer, MMMM As Integer
Dim pivot As Double, temp As Double, arr1() As Variant, arr2() As Variant, x As Integer, y As Integer
Dim determine As Integer, xx As Integer, yy As Integer, cntr As Integer, fudgefactor As Double
determine = 0
x = UBound(AugMatrix, 1): xx = UBound(TempMatrix, 1) ' =Number of ROWS
y = UBound(AugMatrix, 2): yy = UBound(TempMatrix, 2) ' =Number of COLUMNS' =Number of COLUMNS
' Debug.Print x, y, xx, yy
ReDim arr1(1 To 1, 1 To y): ReDim arr2(1 To 1, 1 To y)

For k = 1 To N
' Locate largest matrix element, use as pivot.
pivot = AugMatrix(k, k): P = k

For L = k + 1 To N ' loop each row

If Abs(AugMatrix(L, k)) < Abs(pivot) Then GoTo EndOfLoop

pivot = AugMatrix(L, k)

P = L
EndOfLoop: Next L ' next row
'Debug.Print pivot
' Swap rows
For j = 1 To N + 1
temp = AugMatrix(k, j)
AugMatrix(k, j) = AugMatrix(P, j)
AugMatrix(P, j) = temp
Next j

' Normalize pivot row
For j = 1 To (N + 1)
'If pivot = 0 Then MsgBox "PIVOT is ZERO"
If pivot = 0 Then Debug.Print "TITS",
TempMatrix(k, j) = AugMatrix(k, j) / pivot
Next j

' Do the Gauss elimination.
For i = 1 To N
If i = k Then GoTo EndOfLoop2
For j = 1 To N + 1
TempMatrix(i, j) = AugMatrix(i, j) - AugMatrix(i, k) * TempMatrix(k, j)
Next j
EndOfLoop2: Next i
For i = 1 To N
For j = 1 To N + 1
AugMatrix(i, j) = TempMatrix(i, j)
Next j
Next i
Debug.Print k, pivot
Next k
End Sub

 

[electricpete]

The excel Solver tool works quite good on non-linear equations.

Let's say you have
F(x1) = a1
F(x2) = a2
F(x3)= a3
....
F(xn) = an
Where xi are vectors of multiple independent variables and ai are constants.

Make x1...xn your inputs.
You have a computed value of ai and a target value of ai.
Compute error term ei based on the difference.
Compute sum of squares of errors.
Use solver to miminize sum of squares of errors by varying the inputs (x1..xn)



=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[Denial]

I fully endorse ElectricPete's sugestion about using the Solver.

The only drawback with Solver is that you need to feed it a starting situation that is not "too far" from the final solution (where the determination of what is "too far" depends on how nonlinear your problem is). Since in most engineering problems we have a pretty good idea of the answer before we begin (?) this requirement is not usually all that burdensome.

 

[IDS]

mwalshe - Have you tried the Excel MINVERSE() function?

Can you post some sample data?




Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[electricpete]

MInverse is good for linear systems, but typically not as useful for nonlinear.

Good point by Denial about importance of starting guess for solver solution. There is a refinement to the prodedure I discussed above that usually works quite well to compensate for poor initial guess: We add a user-defined "weight" for each error. Then the quantity to minimize is sum of weighted sum of squares. Proceed as follows:
1 - Initially set all weights to 1.0.
2 - Run solver and accept the solution (it will form initial guess for next rounZ).
3 - Set the weight to 10 on the item that has the highest error in step 2 and set all other weights back to 1.0.
4 - Go back to step 2.

This provides an alternate way to move around in the solution space without directly entering solution guesses (to me it's a lot easier to fiddle with weights than fiddle with solution guesses). If there is an exact solution it, this method will usually find it. Often here is no exact solution (just a best fit) and will find there may be two errors tugging back and forth - set both weights to 10 simultaneously.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

Another refinement is to work with weighted fractional errors instead of just weighted errors. That way the large magnitude resultsw (ai) don't end up have more "weight" than the small ones.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.