sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.84k stars 4.4k forks source link

Singularity Functions that never return infinity for mechanics application #26696

Open mvg2002 opened 3 months ago

mvg2002 commented 3 months ago

The definition of singularity functions of negative order, $< x - a > ^n$ with $n < 0$, states that the function can be evaluated as 0 for all $x$, except when $x = a$, then it is evaluated as infinity. Singularity functions that are evaluated as infinity can cause problems in mechanics calculations.

Right now these singularity functions of negative order are removed from some calculations in the Beam module. An example of how this is done right now:

#Removes the singularity functions of order < 0 from the bending moment equation used in this method
non_singular_bending_moment = sum(arg for arg in self.bending_moment().args if not arg.args[1].args[2] < 0)

This way of handling negative order singularity functions works, but it would be nice if there were special singularity functions for a use in mechanics that make this kind of removing the functions obsolete.

These special mechanics singularity functions would work exacly the same as the normal singularity functions, except that they always return 0 when they are of negative order. In mechanics calculation the singular point of infinity has no meaning and should actually be evaluated as 0.

One of the ways this could be done is by adding a helper method _mechanics_singularity_function in the beam module. This would make the function available for calculations in this module.

An even better way (in my opinion) would be to make the mechanics_singularity_function available to all users. This way users that do mechanics calculations 'by hand' instead of the beam module can also make use of the special singularity functions.

moorepants commented 3 months ago

I just read this independent of the other discussions and I don't think there is enough information here to act on this. More info on the mathematical definition of a singularity function is needed and how the SymPy function does or doesn't adhere to that definition would help. You imply there are two types of singularity functions: what we have and a special one. It isn't so clear what is mathematically meant by special in this context.

mvg2002 commented 3 months ago

First of all I want to state that I do not want to change anything to the current implementation of singularity functions. The current implementation follows the mathematical definition, which it should.

The singularity functions are used a lot in mechanics applications. In these mechanics applications the fact that negative order singularity functions can return infinity can be a problem. This infinity has no meaning in mechanics and can mess up results.

My proposal is to add an extra kind of singularity function that does not return infinity. This can be called something like MechanicsSingularityFunction. This function can then be used in mechanics applications. This function doesn't completely follow the defined rules for singularity functions, but will be much better suited for a use in mechanics applications. This should be clearly stated in the documentation of the function.

To elaborate, the only difference between the 'normal' singularity functions and the one that would be nice for a mechanics application is what they return when they are negative order and $x = a$ in $< x-a > ^n$ with $n<0$.

The normal singularity function returns $\infty$ when $x = a$ in $< x-a > ^n$ with $n<0$. As it should.

The mechanics singularity function would return 0 when $x = a$ in $< x-a > ^n$ with $n<0$.

moorepants commented 3 months ago

Is this connected to the debate of whether a Heaveside evaluated at 0 is 0, 1 or all values in between? Or that the DiracDelta is not infinity at 0?

You second comment block mostly said what you said in the first comment block, "This infinity has no meaning in mechanics and can mess up results." which is not so clear what you mean by "mess up results".

mvg2002 commented 3 months ago

This has nothing to do with the Heaveside function. It does have to do with the DiracDelta function. For mechanics applications the fact that these function can become infinity at a particular point is not wanted.

An example of this:

Let's say this is the bending moment equation for a beam with a hinge: $-EI * P_5 < x-5 > ^ {-1} + F < x-5 > ^1 - R0 < x > ^1 - R{10} < x-10 > ^1 - M_{10} < x-10 > ^{0}$

The hinge is located at $x = 5$. We can use the fact that the bending moment is 0 at the location of the hinge to set up an equation. We fill in bending moment at $x = 5$ = 0. $-EI P_5 \infty - R_0 * 5$.

This is not wanted. The hinge function should not be present in this equation. The $< x-5 > ^ {-1}$ should have become 0 instead of $\infty$.

This problem occurs for quite some calculations in the Beam module. Right now, to sidestep this problem, new equations are made without the negative order singularity function before filling in the values. This works but it would be nicer if the functions could stay in the equation but would just return 0 instead of $\infty$.

moorepants commented 3 months ago

So $< x - a >^{-1}$ is the Dirac Delta function and working with them directly is problematic. Do we need to work with limits to ensure terms behave well when x -> a? My hunch is that just making it zero is not the correct solution and that there must be a mathematically sound way of handling this Dirac Delta.

mvg2002 commented 3 months ago

The problem occurs when substituting values for $x$ in the equations. For example the bending moment equation I gave.

For a hinge at a location $a$ the singularity function that is added to the load equation looks like: $EI * P_a < x - a > ^{-3}$.

After integrating twice this becomes $EI * P_a < x - a > ^{-1}$ in the bending moment equation.

For each hinge the boundary condition is that the bending moment is 0 at the location of the hinge. So at $x = a$.

It will thus always be that we fill in exacly the location of the singularity: $EI P_a < a - a > ^{-1}$ = $EI P_a \infty$, which should be $EI P_a * 0$ = 0.

I must say that my knowledge on working with limits is not great, but I can't immediatly think of a way they can help here. The problem is that the location where the singularity funcitons is evaluated is always exactly the location of the singularity.

moorepants commented 3 months ago

Does "which should be $EI P_a 0$ = 0" mean "if it were this, then my problem is solved"? "Should" implies it should be that mathematically. That I'm not sure of.

Maybe the limit needs to be applied before integrating or something. I'd have to review the mathematics more carefully to give any useful answer.

mvg2002 commented 3 months ago

Maybe I stated it to harsch. What I mean is $EI P_a < x - a > ^{-1}$ should no more be present in the equation after applying $x=a$.

In the example I gave:

After filling in $x=a$, the equation should be $-R_0 5$ instead of $-EI P_a \infty - R_0 5$