py_UMom 1.04

The py_UMom spreadsheet has been updated with further revisions to the py_OptShearCap3600 function. The new version can be downloaded from:

py_UMom.zip

For more information on the other functions in the spreadsheet see: py_UMom spreadsheet and OptShearCap3600 function.

The changes to the py_OptShearCap 3600 function are:

  • There is a new option to calculate the longitudinal force due to shear using the AS 5100.5 equation, with all other calculations to AS 3600.
  • The input has been modified so that axial force, moment, shear and torsion are specified for the critical load case, and a range of M/V (moment/shear ratio) values are entered to generate the interaction diagram.
  • The option to adjust axial load has been removed. The input axial load will be applied to all load cases.

The revised input and typical output is shown in the screenshots below:

Examples comparing the different longitudinal force options with different levels of shear reinforcement are shown below. The 4 options in each graph are:

  • Longitudinal shear force to AS 3600 with no adjustment of the compression strut angle, Theta.
  • As above with the compression strut adjusted to minimise the longitudinal shear force.
  • Longitudinal shear force to AS 5100.5 with the vertical force in the shear steel, Vus, limited to V*.
  • As above with no limit to Vus.

The first case has shear reinforcement just adequate for the input V*:

Increasing the shear reinforcement greatly increases the shear range with unconservative unreduced bending capacity when the AS 5100 calculation is used with no limit on Vus:

When the shear reinforcement is reduced below the level required for the input shear force the critical forces are all in the range where shear capacity governs the failure mode, and all four options give the same shear capacity:

Posted in Beam Bending, Concrete, Excel, Link to Python, Newton, PyXLL, UDFs | Tagged , , , , , , , | 1 Comment

Mandelbrot, Mojo, Numpy and Numba

I recently read an article about connecting Python to the new Mojo language:

Python Can Now Call Mojo

The article had some code examples, including code to plot the Mandelbrot set, comparing plain Python, Python with Numpy, and Mojo, and what particularly interested me was that the code with Numpy was about 4 times faster than plain Python, whereas I had found Numpy for this task was slower than Python, as reported here: Computing and plotting the Mandelbrot set in Excel …

I have now updated the matplot-mandelbrot spreadsheet to include revised Numpy code, based upon on the code at the link above, and found that the Numpy code went from up to about 2 times slower to 3-4 times faster.

The updated spreadsheet and Python code can be downloaded from: FastMandelbrot.zip

The Numba code using Guvectorize remains by far the fastest, but the revised Python + Numpy code is now much faster than the plain Python. The original code loops through the results array cell by cell, converting the coordinates to a complex number then calling a plain Python routine to do the calculations:

def mandelbrot_set1(xmin,xmax,ymin,ymax,width,height,maxiter):
    r1 = np.linspace(xmin, xmax, width)
    r2 = np.linspace(ymin, ymax, height)
    n3 = np.empty((width,height))
    for i in range(width):
        for j in range(height):
            n3[i,j] = mandelbrot(r1[i] + 1j*r2[j],maxiter)
    res = (r1,r2,n3)
    return res

The revised code generates an array of complex numbers in a single operation, then generates the Mandelbrot values with a Numpy vectorized computation:

def mandelbrot8(x, y, maxiter, out = 0):
    """Generates a Mandelbrot set using NumPy for vectorized computation."""
    c = x[:, np.newaxis] + 1j * y[np.newaxis, :]
    z = np.zeros_like(c, dtype=np.complex128)
    image = np.zeros(c.shape, dtype=int)
    
    for n in range(maxiter):
        not_diverged = np.abs(z) <= 2
        image[not_diverged] = n
        z[not_diverged] = z[not_diverged]**2 + c[not_diverged]
        
    image[np.abs(z) <= 2] = maxiter
    return image  

The spreadsheet plots the calculated shape with chosen colour properties to full scale:

or scaled up by 100,000:

or 1,000,000:

or by 10,000,000:

or by 100,000,000:

and more, but after 1 billion, the image starts to lose resolution using 64 bit floats.

For those who want to go further, visit YouTube:

“The deepest Mandelbrot zoom ever! This zoom is over 24-hours long, uploaded to YouTube in 11 parts. It will finish on a mini-Mandelbrot. As such, I am claiming a world record for the deepest Mandelbrot video! I have kept to the traditional colouring style for this one. The first video may seem a little repetitive at first, but this builds more interesting shapes later into the zoom. You’re welcome to skip ahead a little. Consider playing your own music while you let this one roll. Sit back, relax, and soak it in. (The last minute also generates a cool optical illusion with the pattern it finishes on).”

Posted in Arrays, Drawing, Excel, Link to Python, Maths, Newton, NumPy and SciPy, PyXLL, UDFs | Tagged , , , , , , , | Leave a comment

py_Numpy update 2: Polynomial functions

Continuing from the previous post, this post looks at the polynomial functions in the revised py_Numpy spreadsheet. The updated code and spreadsheet can be downloaded from:

 py_SciPy.zip

The available functions call the Numpy functions for solution and evaluation of polynomials, plus some additional functions with closed-form solutions for polynomials of up to quartic order. The available functions for finding the roots of polynomial equations are shown below (click on any image for full-size view):

Numpy iterative functions:

  • py_PolyRoots: Find the real and complex roots of a polynomial equation, calling the Numpy polyroots function. Coefficients are listed in increasing order of x.
  • py_PolyRoots: Provides an option for the order of coefficients. By default list in decrasing powers of x
  • py_PolyRootsC: Accepts complex values for the input coefficients. Also has an option to list coefficients in descending powers of x, default = ascending powers.
  • py_Polyroots2: Find the real and complex roots of a polynomial equation, calling the Numpy Polynomial function and .roots method.

Closed-form functions:

  • py_Quadratic: Find the real and complex roots of a quadratic equation: a x^2 + b x + c = 0
  • py_Cubic: Find the real roots of a cubic equation: a x^3 + b x^2 + c x + d = 0
  • py_CubicC: Find the real and complex roots of a cubic equation: a x^3 + b x^2 + c x + d = 0
  • py_CubicT: Find the real roots of a cubic equation: a x^3 + b x^2 + c x + d = 0, alternative solution method.
  • py_Quartic:  Find the real and complex roots of a quartic equation: a x^4 + b x^3 + c x62 + d X + e = 0

Calculation times for 100,000 iterations are shown in the screenshot below, with plain Python code, and with the Numba JIT compiler.

More functions:

  • py_PolyFromRoots: Return the polynomial with the given roots
  • py_PolyVal Evaluate a polynomial for a given x value
  • py_PolyFit Fit a polynomial function to XY data
  • py_PolyCompanion Returns the companion matrix to an array of polynomial coefficients in ascending order
  • py_PolyDer Differentiate a polynomial
  • py_PolyInt Integrate a polynomial
  • py_PolyAdd Add one polynomial to another
  • py_PolySub Subtract one polynomial from another
  • py_PolyMul Multiply one polynomial by another
  • py_PolyDiv Divide one polynomial by another
  • py_PolyPow Raise a polynomial to a power
Posted in Curve fitting, Excel, Link to Python, Maths, Newton, NumPy and SciPy, PyXLL, UDFs | Tagged , , , , , , | Leave a comment

Numpy where and split, and py_Numpy update

The py_Numpy spreadsheet and associated code has been updated, with two new functions, links to the Numpy on-line help updated, and updates to sorting and polynomial functions. The updated code and spreadsheet can be downloaded from:

Updated Numpy files included in: py_SciPy.zip

Note that the Numpy code comes in 2 versions: pyNumpy-noJIT.py and pyNumpy-jit.py. The pyNumpy-jit.py requires the Numba just-in-time compiler to be installed, which provides very much faster results for maths intensive operations, compared with plain Python code. The default pyNumpy.py module is currently a copy of the noJIT version.

The py_SciPy download contains full open-source Python code and example spreadsheets for a wide range of Scipy functions, as well as the Numpy functions covered in this post. For more details, see the series of posts starting at: Scipy functions with Excel and pyxll. Also see Python and pyxll for details of the pyxll add-in, required to link my Python code to Excel.

This post will look at the two new functions, py_Where and py_Split, which can be found on the Array Functions tab of the py_Numpy spreadsheet. The other changes will be covered in the following post.

The py_Where function calls the Numpy “where”:
numpy.where(condition, [xy, ]/)
Return elements chosen from x or y depending on condition.

The x, y return values are optional, and where omitted the function returns an array with the base 0 index values of the data values that satisfy the given condition, as shown in Column D below:

For this case the required inputs are:

  • the data range (a single column or row)
  • the condition, entered as text, one of <, <=, =, >=, or >
  • the limit value

The output is a single column array with the base 0 index of the data values that satisfy the condition, in this case those that are <= 13.

Where x and y are specified typical output is shown in Column G above. Where the condition is satisfied x (‘Hearts’) is returned, and otherwise y (‘Spades’).

The py_Split function splits the data into a number of smaller arrays at the specified locations:

The sub-arrays are returned as rows, with empty cells displaying #N/A.

The Numpy help on the function (and most of the other functions in py_Numpy) can be displayed by clicking on ‘Help on this function’ in the bottom-left corner of the function dialog box:

Posted in Arrays, Excel, Link to dll, Link to Python, PyXLL, UDFs | Tagged , , , , , , | 1 Comment

py_RC Elastic 1.05 with elastic biaxial bending

The py_RC Elastic spreadsheet (last discussed here) has been updated with a new function py_Biax, providing elastic analysis of reinforced concrete sections under combined axial load and biaxial bending. The new spreadsheet and associated Python code can be downloaded from:

py_RC Elastic.zip

The download file includes full open-source Python code. For details of the pyxll package required to link the Python code to Excel see: https://newtonexcelbach.com/python-and-pyxll/

The code and spreadsheet layout are based on the VBA version presented at: Elastic Biaxial Bending

Input for a T-section beam is shown in the screenshot below. The concrete is defined by the coordinates of each corner point, and reinforcement is is detailed in layers, with the coordinates of the ends of each layer:

The angle of the neutral axis is found by iteration, and the estimated angle should be entered in cell E5, then click the “Adjust NA Angle” button. The NA angle and position will be adjusted so that the stress at both ends of the NA is zero:

The function output displays:
Concrete, reinforcement and combined section properties:

Concrete stresses at each corner of the section in compression:

Reinforcement stresses at the ends of each layer:

Posted in Beam Bending, Concrete, Excel, Link to Python, Newton, PyXLL, UDFs | Tagged , , , , , , | Leave a comment