Solving non-linear equations with two or more unknowns – 5

This will be the last of the series on solving non-linear equations (for now).  Up until now all the examples have had two unknown values, and two target values.  This can be extended by making three changes to the code:

  • Set up an nxn matrix of the function slopes with respect to each of the unknown values (the function Jacobian)
  • Calculate the first estimate of the target values, using an estimated value for each of the unknowns.
  • Form and solve a series of simultaneous equations to find a better estimate of the unknowns.

In Excel VBA the solution to the equations can be found using the Worksheetfunction.MInverse and MMult functions, as shown in the code below:

LoopNum = LoopNum + 1
        ' Evaluate function at guessed values and small increments.
        For i = 1 To NumVar
            Var2A(i, 1) = Var1A(i, 1) * VarFact
            VarDiffA(i, 1) = Var1A(i, 1) * (VarFact - 1)
        Next i
        Res1 = Application.Run(Func, Var1A, X)
        For j = 1 To NumVar
            ResDiff(j, 1) = Target(j, 1) - Res1(j, 1)
        Next j
        Temp = Var1A(1, 1)
        Var1A(1, 1) = Var2A(1, 1)
        For i = 1 To NumVar
            Res2 = Application.Run(Func, Var1A, X)
            For j = 1 To NumVar
                SlopeA(j, i) = (Res2(j, 1) - Res1(j, 1)) / VarDiffA(i, 1)
            Next j
            Var1A(i, 1) = Temp
            If i < NumVar Then
                Temp = Var1A(i + 1, 1)
                Var1A(i + 1, 1) = Var2A(i + 1, 1)
            End If
        Next i
        ' Solve SlopeA
        InvA = WorksheetFunction.MInverse(SlopeA)
        ResA = WorksheetFunction.MMult(InvA, ResDiff)
        ErrSum = 0
        For i = 1 To NumVar
            ErrA(i, 1) = (Target(i, 1) - Res1(i, 1))
            ErrSum = ErrSum + Abs(ErrA(i, 1))
            Var1A(i, 1) = Var1A(i, 1) + ResA(i, 1)
        Next i

The User Defined Function (UDF) MSolve has been modified as shown above, and the earlier version renamed MSolve2.  A spreadsheet including examples of both functions, and full open source code, can be downloaded from:

Input and output from both functions is shown in the screenshot below:


MSolve2 works as before, in the example returning the depth of neutral axis and compression face strain in a reinforced concrete section, for any given axial load and bending moment.

The MSolve example adds one more unknown (the diameter of the tension reinforcement), and another target, the tension reinforcement stress.

The Python Scipy package also contains a number of routines for solving problems of this type.  In order to access these from Excel we need:

  • A python routine for the function to be solved.
  • An interface allowing this routine and the Scipy functions to be accessed from Excel.

I have updated the Eval-PyIntegration spreadsheet to perform this task, and this can be downloaded from:

The spreadsheet uses ExcelPython for communication with Python, and the download file includes all necessary files (other than Excel, Python and Scipy).

The example below shows the xl_SolveF function used to find the depth of neutral axis, compression face strain, and tension bar diameter required for a specified axial force, bending moment and tension steel stress.  In this example the Python function being solved (RCForceMS) is the simplified version, ignoring concrete in tension.


The download also includes Python versions of the functions CurveatMA and CurveatMAS, with input as shown below.  The functions make use of the SciPy optimize.root function, which has many optional arguments, detailed in the SciPy manual.  In the spreadsheet the only options included are for the solver type, and tolerance for termination.


The screenshot below shows the results of CurveatMA for an axial load of zero and bending moment of zero.  The returned values of depth of neutral axis (DNA) and compression face strain (Epsc) have been checked using the function SectForceMV, which finds that the axial force and bending moment are as expected.

The xl_SolveF function can also be called directly, using a function that returns the difference of the function values from the target values, in this case SectForceMDiff.


Finally, CurveatMAS adjusts the neutral axis depth, compression face strain, and tension steel area such that the axial force, bending moment and section curvature match the target values.  Once again, xl_SolveF may be used with SectForceMdiff to return the same 3 values.

Note that the steel force and tensile strain found by these functions are average values, including tension stiffening effects.  The results are intended for use in deflection calculations, and should not be used in strength or stress limit calculations.




This entry was posted in Beam Bending, Concrete, Excel, Link to Python, Maths, Newton, NumPy and SciPy, UDFs, VBA and tagged , , , , , , , , , , , , . Bookmark the permalink.

3 Responses to Solving non-linear equations with two or more unknowns – 5

  1. lori says:

    Doug – nice series of posts… i have found convergence is very fast with this method (around 10 digits of precision for only 5 iterations with only fairly approximate initial values).

    A possible tweak is to substitute .MMult(.MInverse(SlopeA), ResDiff) for .Linest(ResDiff, SlopeA, 0) which returns equivalent results but in reverse order. Linest has become more accurate since version 2010 and no longer has limitations on matrix size, it also allows for near singular or non-square matrices which then makes this approach equivalent to the Gauss-Newton method.

    I saw Georg referred to Levenberg-Marquardt in an earlier comment which is a modification for when convergence may be an issue although there’s a trade off as it may end up being slightly slower.


    • dougaj4 says:

      Thanks Lori, very interesting. I’ll have a look at incorporating the Linest option when I have time.

      Regarding the Levenberg-Marquardt method, that’s one or the options of the SciPy function. When I first started with that I was passing a depth value in mm (of the order of 100) and a strain value of the order of 1E-4, and the Levenberg-Marquardt (set Solver Type to “lm”) was the only one that worked. There are various options for setting limits and tolerances which probably also would have worked with the other solvers, but I just changed the strain input to microstrain, so it was also of the order of 100, and with that change the hybrid solver also converged well, and was quite a bit faster.

      For my current application the Python/Scipy solvers don’t offer any real advantage over my VBA (they are actually quite a bit slower), but I am sure there are many applications where the additional solver options would be useful.


  2. Pingback: xlSciPy – Python SciPy for Excel; Update with new functions | Newton Excel Bach, not (just) an Excel Blog

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s