Continuum Mechanics's Beam - convert SingularityFunction to Piecewise for plotting?

[Davide-sd]

The sympy.physics.continuum_mechanics.beam uses the SingularityFunction to create symbolic expression for the loads. The plotting module converts the symbolic expression to a lambda function to be evaluated with NumPy. However, Numpy has no idea what SingularityFunction is, hence it fails. Then, the plotting module evaluates the expression with SymPy, which is slow.

Converting the symbolic expression to a Piecewise allows the plotting module to use NumPy (potentially better performance and better quality). I wonder if it's always possible to convert a SingularityFunction to a Piecewise in the context of the continuum mechanics module.

Here is a little benchmark:

from sympy import *
x = symbols("x")
expr1 = 13750*SingularityFunction(x, 0, 0) - 5000*SingularityFunction(x, 2, 0) - 10000*SingularityFunction(x, 4, 1) + 31250*SingularityFunction(x, 8, 0) + 10000*SingularityFunction(x, 8, 1)
# same expression but containing Piecewise
expr2 = piecewise_fold(expr1.rewrite(Piecewise))

I'm going to generate numerical data for plotting using adaptive=True and adaptive=False, for both expressions. For an ordinary Matplotlib plot, the width in pixels covered by a curve is between 500-600px. For the uniform meshing strategy, I use n=1000 discretization points. It means there will be (on average) 1000points/550px = 1.818181 points/px in the horizontal direction, which is more than enough to properly capture discontinuities (jumps in the y-direction).

from sympy.plotting.series import LineOver1DRangeSeries
# evaluates SingularityFunction with SymPy
s1a = LineOver1DRangeSeries(expr1, (x, 0, 8), adaptive=True)
s1b = LineOver1DRangeSeries(expr1, (x, 0, 8), adaptive=False, n=1000)
# evaluate Piecewise with NumPy
s2a = LineOver1DRangeSeries(expr2, (x, 0, 8), adaptive=True)
s2b = LineOver1DRangeSeries(expr2, (x, 0, 8), adaptive=False, n=1000)

These are the results (computed with a modern high end cpu):

# evaluates SingularityFunction with SymPy
%timeit -r 5 -n 50 s1a.get_data()
# 49.5 ms ± 110 µs per loop (mean ± std. dev. of 5 runs, 50 loops each)
%timeit -r 5 -n 50 s1b.get_data()
# 272 ms ± 390 µs per loop (mean ± std. dev. of 5 runs, 50 loops each)

# evaluate Piecewise with NumPy
%timeit -r 5 -n 50 s2a.get_data()
# 12 ms ± 137 µs per loop (mean ± std. dev. of 5 runs, 50 loops each)
%timeit -r 5 -n 50 s2b.get_data()
# 56.5 ms ± 1.03 ms per loop (mean ± std. dev. of 5 runs, 50 loops each)

The last and first timings show that it would be possible to plot it with adaptive=False and Numpy at roughly the same "speed" as adaptive=True with SymPy, but potentially increasing the quality of the plot. The second timing shows that evaluating an expression containing SingularityFunction with the uniform meshing strategy is out of the question: it would take well over one second to generate a plot with all the loads.

For example, a plot of an expression containing SingularityFunction with adaptive=True evaluated with SymPy:

plot(expr1, (x, 0, 8), adaptive=True)

There is a discontinuity at x=4 and it would appear there is another discontinuity at roughly x=4.2: that's an artifact of the adaptive algorithm, which has reached its maximum depth of recursion. Note that by running this command again, the artifact could disappear because the adaptive algorithm uses random numbers. One could be tempted to use a seed in the adaptive algorithm in order to always get the same results. However, that's not going to resolve the depth of recursion limit. Using adaptive=False eliminates the artifacts, producing a better plot taking the same computational time:

plot(expr2, (x, 0, 8), adaptive=False, n=1000)

No artifacts at x=4.2.

[oscarbenjamin]

However, Numpy has no idea what SingularityFunction is, hence it fails. Then, the plotting module evaluates the expression with SymPy, which is slow.

Maybe lambdify should convert to Piecewise.

[moorepants]

The singularity functions should always be able to be rewritten as piecewise. These functions are a shorthand for piecewise functions, in some sense.

[Davide-sd]

Maybe lambdify should convert to Piecewise.

This would be awesome.