A Hard Rain

I happened to be looking at an old post here; Various routes to and from Nottamun Town, featuring Bob Dylan’s performance of A Hard Rain’s Gonna Fall, which led me to Patti Smith’s performance of the same song at last year’s Nobel Prize ceremony, and this is a unique performance in very many ways:

Advertisements
Posted in Bach | Tagged , , , | Leave a comment

Using Flash Fill

There have been several posts here on using functions or user defined functions (UDFs) to extract text or numbers from a longer string (most recently, here), but for many purposes it is quicker and easier to use the “flash-fill” feature, introduced in Excel 2013.

The examples below show examples of extracting a number from the end of a text string, using data copied and pasted from Strand7 as an example:

In the first example it is required to extract the numbers from the text in column B.  Simply enter the first number in the adjacent column, then press Ctrl-E (or click Data-Flash Fill):

The column is immediately filled with the required values.  A small icon appears next to the top cell, which gives the options displayed above if selected.

It is also possible to use a non-adjacent column:

Note that in this case the values extracted are not in ascending order.

More complex examples also work:

In this case the values in rows 14 to 17 (5-1 to 5-4) have been extracted as a date.  This can be fixed by simply formatting the example cell as text:

Also if the example cell is changed immediately, the flash-fill automatically updates:

There are however times when more control is required.  The screen-shot below shows results from flash-fill (in column A), compared with the ExtractNums UDF in Column D (download here).

Posted in Excel, UDFs, VBA | Tagged , , , | 1 Comment

Hofstadter, Pahud, Bach

The latest ACO concert series features the Swiss flautist Emmanuel Pahud (interview here). The performance included J.S. Bach’s 6 Part Ricercar, from The Musical Offering, which was the inspiration for the final dialogue in Douglas Hofstadter’s Godel, Esher, Bach:

Achilles has brought his cello to the Crab’s residence, to engage in an evening of chamber music with the Crab and Tortoise. He has been shown into the music room by his host the Crab, who is momentarily absent, having gone to meet their mutual  friend the Tortoise at the door. The room is filled with all sorts of electronic equipment-phonographs in various states of array and disarray, television screens attached to typewriters, and other quite improbable-looking pieces of apparatus.  Nestled amongst all this high-powered gadgetry sits a humble radio. Since the radio is the only thing in the room which Achilles knows how to use, he walks over to it, and, a little furtively, flicks the dial and finds he has tuned into a panel discussion by six learned scholars on free will and determinism. He listens briefly and then, a little scornfully, flicks it off. …

On the text:

This Dialogue is an exuberant game played with many of the ideas which have permeated the book. It is a reenactment of the story of the Musical Offering, which began the book; it is simultaneously a “translation” into words of the most complex piece in the Musical Offering: the Six-Part Ricercar. This duality imbues the Dialogue with more levels of meaning than any other in the book. Frederick the Great is replaced by the Crab, pianos by computers, and so on. Many surprises arise. The Dialogue’s content concerns problems of mind, consciousness, free will, Artificial Intelligence, the Turing test, and so forth, which have been introduced earlier.  It concludes with an implicit reference to the beginning of the book, thus making the book into one big self-referential loop, symbolizing at once Bach’s music, Escher’s drawings, and Gödel’s Theorem.

From: Gödel, Escher, Bach Wiki 

And the music:

 

Posted in Bach, Maths, Newton | Tagged , , , , , | Leave a comment

Working with implied units – update

Following comments at Eng-Tips I have added two examples with different procedures for evaluating ACI equations for the flexural tensile strength of concrete here.

Posted in Excel, Maths, Newton, UDFs, VBA | Tagged , , , , , , , , , , , , | Leave a comment

Using Alglib least-squares solvers

Following the last post on using the Python version of Alglib from Excel, via xlwings, this post looks in more detail at alternatives for fitting a non-linear function to a set of data, using the Levenberg-Marquardt method.  The spreadsheet with examples and full open-source code may be downloaded from:

xlAlglib.zip

The spreadsheet requires xlwings, Python (including Numpy and Scipy) and Alglib to be installed.  All are available for free download, see Excel to Alglib via xlwings for more details.

In the examples below a non-linear function has been fitted to data generated using the same function plus some random scatter using two different Alglib modules:

  • The lsfit module (least-squares fitting)
  • The minlm module (Levenberg-Marquardt optimizer)

For each of these modules there is one user defined function (UDF) that uses a function entered as text on the spreadsheet, and another that calls a python function.

The lsfit versions (  xl_NLFitwFunc and  xl_NLFitwText) are described with examples on the Fit NonLin sheet:

The minlm versions (xl_minlm_vb and xl_minlm_vbtext) are found on the LM Solvers sheet:

The example shown below uses all four functions to fit a function of the form:
a * exp(b*x) + c * x + d
to the data shown in the screenshot below:

The function xl_minlm_vb (shown below) has the following arguments:

  • x (B42:B45), Initial guess of function parameters, single column range
  • xydat (B22:C34), Data to be fitted, two column range
  • Func (C37), Fitting function; a function in NLFuncs.py
  • bndl (B48:E48), Lower bound values, single row range
  • bndu (B49:E49), Upper bound values, single row range
  • epsx (B51), Maximum error
  • maxits (B52), Maximum iterations
  • out (0), Out = 0 returns results, anything else returns termination type code

Results of the four functions are shown below:

Code for the xl_minlm_vb function is described below.  Open source code for all the others is available in the download zip file.

The function xl_minlm_vb (automatically generated by xl wings) calls the Python function of the same name:

@xw.func
@xw.arg('x', ndim=1)
@xw.arg('bndl', ndim=1)
@xw.arg('bndu', ndim=1)
@xw.ret(transpose= True)
def xl_minlm_vb(x, xydat, funcname, bndl, bndu, epsx, maxits, out = 0):
    """
        Minimise a named function with Alglib minlmcreatev  Levenberg–Marquardt algorithm
        Initial guess of function parameters, single column range
        Data to be fitted, two column range
        Fitting function; a function in NLFuncs.py
        Lower bound values, single row range
        Upper bound values, single row range
        Maximum error, 0 = automatic
        Maximum iterations, 0 = unlimited
        out = 0 returns results, anything else returns termination type code
    """
    if type(x) == tuple: x = list(x)
    if type(bndl) == tuple: bndl = list(bndl)
    if type(bndu) == tuple: bndu = list(bndu)
    nx = len(x)
    nf = len(xydat)
    func = getattr(NLFuncs, funcname)
    state = xal.minlmcreatev(nx, nf, x, 0.00001)
    xal.minlmsetbc(state, bndl, bndu)
    xal.minlmsetcond(state, epsx, maxits)
    xal.minlmoptimize_v(state, func, None, xydat)
    x, rep = xal.minlmresults(state)
    if out == 0:
        return x 
    else:
        return rep.terminationtype
  • if type(x) == tuple:  … Arrays are passed from Excel to Python as tuples, but Alglib requires a list.
  • func = getattr(NLFuncs, funcname): looks in NLFuncs.py for a Python function with the name given by funcname, and assigns this function to “func”.
  • state = xal.minlmcreatev( … x, rep = xal.minlmresults(state):  calls Alglib minlm functions to adjust the parameters passed to func in the “state”  object to minimize the square of the results.

“func” in this case is xl_expfunc:

def xl_expfunc(x, fi, func_dat):
    nf = len(func_dat)
    func_dat  = np.array(func_dat) 
    for i in range(0,nf):
        X = func_dat[i,0]
        fi[i] = x[0]*np.exp(x[1]*X)+x[2]*X+x[3] - func_dat[i,1]
    return

Note that Alglib minimises the square of the array returned by xl_expfunc by adjusting the values in the x array.  To use this for curve fitting purposes we have to subtract the supplied Y value from the exponential function:
fi[i] = x[0]*np.exp(x[1]*X)+x[2]*X+x[3] – func_dat[i,1]

Alglib adjusts the function parameters until the error or maximum iterations criteria are satisfied, then extracts the required results from the “state” object with:
x, rep = xal.minlmresults(state)

The x array is then returned to the spreadsheet, via the VBA function.

Posted in AlgLib, Curve fitting, Excel, Link to Python, Maths, Newton, NumPy and SciPy, UDFs, VBA, xlwings | Tagged , , , , , , , , , , | Leave a comment

Working with implied units

The Units4Excel spreadsheet provides a number of functions to convert between different units and evaluate functions entered as text on the spreadsheet, taking account of the units of both the inputs and the results.  The latest version can be downloaded from:

Units4Excel.zip

I have now added some examples comparing reinforced shear capacity to the ACI 318 design code, using the “customary US units” and SI units versions of the code:

The examples shown below use the user defined functions (UDFs) Eval and EvalU.  Eval returns the value of any function entered as text on the spreadsheet, substituting given values for any parameters.  EvalU is a “unit aware” version, returning the result in base SI units, or in any specified units.  EvalU returns the calculated units if they are not specified, or an error message if the specified output units are not consistent with those calculated from the input units.

Evaluating formula 11-4 as given in the code produces strange results, even if consistent input units are used.  If the input units are ignored (using the Eval function) the formula returns the correct results, provided that the input has consistent units, but using EvalU with no specified return unit finds a unit with non-integer exponents, and including m, rather than just N.  If the return units are specified, as a force unit, this results in an error message, since the specified return units are different to those calculated from the input units:

The problem is that the code formula contains implied units, since Nu/(2000*Ag) has units of stress (rather than no units), and fc^0.5 has units of stress^0.5, rather than stress.  The units can be made consistent by converting the 2000 (or 14 for the SI version) to an input parameter, Kc1, with stress units, and adding a new parameter, Kc2, also with stress units, as shown below:

EvalU now returns calculated units of N, as expected, and both versions of the formula return very similar results in N.  Also note that when the output units are consistent with the input units (either lbf or N), EvalU and Eval return exactly the same value, as would be expected.

It is also possible to incorporate the initial constants (2 or 0.17) in the Kc1 and KC2 factors, and to use inconsistent input units.  In the first example below Kc1 and Kc2 have been changed to ksi units, and Kc2 is factored by 2^2, so the initial 2 factor can be removed.  In the SI example the units of Kc1 and Kc2 are unchanged, but Kc2 is factored by 0.17^2:

With the revised formulas the EvalA results are correct, but the Eval results, ignoring the input units, are now incorrect, because the input units are inconsistent (stresses are in ksi or MPa, but forces are in lbf or kN).

Formula 11-6 has no factors with implied units, so can be evaluated with no additional parameters, but as before, the input and output units must be consistent if Eval is used.  EvalU will return correct results with any recognised moment output units:

The results of 11-6 can now be used to evaluate 11-5 for sections with combined axial load and shear:

Finally the upper limit to the shear force for combined axial load and shear is found with formula 11-7.  Note that in this case the final expressions in the formula have a different form in the two versions of the code (1 + Nu/(500*Ag) for the US units version and 1 + 0.29*Nu/Ag for the SI version).  Both have been converted to the form 1 + Nu/(Kc2*Ag), so Kc2 can have units of stress, so for the SI version the value of Kc2 becomes 1/0.29:

Two additional points to note when using EvalU:

  • The units are case sensitive, so if Excel helpfully insists on correcting MPa to Mpa, the units will not be recognised.  To avoid this problem go to File-Options-Proofing, click AutoCorrect Options, and either unselect “Correct TWo Initial CApitals, or click Exceptions and add MPa to the exceptions list on the Initial CAps tab.  For a quick one-off you can also enter a space after MPa, then press Ctrl-Z.
  • Both input units and output units (if used) must be recognised units with the correct abbreviation.  The ListNonSI UDF will list all recognised units for any specified unit type, for example “area” as shown in the screen-shot below.  To display all the results enter the function with the unit type, then press Ctrl-Shift-S to resize the output array.

Update 30 Sep 2017:
Following comments at Eng-Tips I have added two examples with different procedures for evaluating ACI equations for the flexural tensile strength of concrete:

The screenshot below shows equations from the US units and SI units versions of ACI 318, evaluated using the procedures described above, i.e. the concrete compressive strength has been factored by one, with units appropriate to the code.  Note that the SI results are about 12% higher than the US units results, because the Kc1 factor is rounded up in the SI code:

An alternative approach is to divide the fc value by the units specified for the equation.  The initial factor should than be given units of stress, or alternatively multiplied by the alternative stress units.  Unlike Mathcad (which recognises units within a formula), EvalU requires the units to be specified as a parameter, with a value of 1, as shown below:
The third example shows the US units formula evaluated with fc in MPa, and return units specified as psi, showing exact agreement with the US units version.

Posted in Concrete, Excel, Maths, Newton, UDFs, VBA | Tagged , , , , , , , , , , , , | 2 Comments

Stepping through 3DFrame

Following a recent comment this post looks at the code for setting up the stiffness matrix for a 3D frame, and how to step through any selected part of the code.

The general procedure for setting up a 3D stiffness matrix, and rotating it to the global coordinate system is given at the link above.  The 3DFrame spreadsheet (including full open-source code) can be downloaded from:

3DFrame.zip

The main program for the 3D analysis is the sub-routine p44, found in the mFrameP44 module in the Visual Basic Editor (VBE), which can be accessed by:

  • From the spreadsheet, press Alt-F11, or click the Visual Basic icon on the Developer tab*
  • Double-click the mFrameP44 module in the Project Explorer window.
  • Click the right-hand dropdown box above the editor window, and select P44

*  If the Developer tab is not visible, go to File-Options-Customize Ribbon, and select it as one of the main tabs.

The code for setting up the global stiffness matrix is:

    ' !---------------------global stiffness matrix assembly------------------------
elements_2:
' Range("kga").ClearContents
    For iel = 1 To Nels
        For i = 1 To nod
            Num(i) = g_num(i, iel)
            For j = 1 To NDim
                Coord(i, j) = g_coord(j, Num(i))
            Next j
        Next i

        ReDim km(ndof, ndof)
        i = BeamConER(iel, 6)
        If i = 0 Then
            ReDim R0(1 To 3, 1 To 3)
            Call rigid_jointed3(km, Prop, Gamma, Etype, iel, Coord, R0)
        Else
            For j = 1 To 12
                SpringA(1, j) = BeamHinge(i, j + 1)
            Next j
            Call spring_jointed3(km, Prop, Gamma, Etype, iel, Coord, R0, SpringA)
        End If

        For i = 1 To ndof
            g(i) = g_g(i, iel)
        Next i

        Call formkv(kv, km, g, neq, iel, g_coord, Num)
    If iel = Range("kmbeam").Value2 Then Range("km") = km
    Next iel    ' elements_2

This code loops through each beam element in the frame and calls either the rigid_jointed3 or spring-jointed3 routines for beams with fixed or spring end conditions respectively. The 12×12 matrix returned for each beam is then added to the global stiffness matrix with the formkv routine.

To step through the code in the p44 routine simply make sure that the cursor is somewhere within the body of the routine and press F8, or Debug-Step Into.  To run the routine to a selected point either:

  • Set a breakpoint at the required location by clicking in the left hand margin, then press F5 or Run-Run or
  • Select the line to run to, then press Ctrl-F8 or Debug-Run to cursor.

There is also a separate function, rigid_jointed3TT, that can be called from the spreadsheet for a specified beam, and returns the beam 12×12 stiffness matrix to the spreadsheet as an array function.  The instructions below show how to step through this function line by line.

To step through a VBA function it must either be called from a sub-routine, or initiated from the spreadsheet.  To do the latter:

  • Create a breakpoint somewhere within the VBA code by clicking in the left-hand margin.
  • Select a cell in the spreadsheet that calls the chosen function.  The screenshots below use cell C35 in the Matrix examples spreadsheet.
  • Press F2 (edit), then either enter for a function that returns a single value, or in this case Ctrl-Shift Enter to return an array.  The function will then run to the selected breakpoint:

You can now step through the code one line at a time by pressing the F8 key (or Debug-Step Into).  Note that the Locals Window, under the main Editor window displays the status and value of all variables, and you can also see the value of a variable by simply placing the cursor over the variable name anywhere in the code:

Procedures for stepping through the code are then the same as for a sub-routine.  To run to a selected line press Ctrl-F8 or right click and Run to cursor:

At the end of the function it will either return the end result (the km2 array), or one of the four intermediate arrays, depending on the value of the Out argument:

The Debug menu provides further options, for example allowing control to jump directly to a specified point, or for called routines to be run in a single step.

For more details of VBA debugging options see Chip Pearson’s Debugging VBA.

Posted in Beam Bending, Excel, Finite Element Analysis, Frame Analysis, Newton, UDFs, VBA | Tagged , , , , , | Leave a comment