## Non linear equations solver using VBA

## Non linear equations solver using VBA

(OP)

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

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

## RE: Non linear equations solver using VBA

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.

## RE: Non linear equations solver using VBA

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.

## RE: Non linear equations solver using VBA

Can you post some sample data?

Doug Jenkins

Interactive Design Services

http://newtonexcelbach.wordpress.com/

## RE: Non linear equations solver using VBA

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.

## RE: Non linear equations solver using VBA

=====================================

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