Repository: sympy/sympy
Pull request: Fixed Maximum bending moment bug in physics beam
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #24221

#### Brief description of what is fixed or changed
Simplification of shear force equation was required before converting it into Piecewise. Also corrected sorting singularity list which was causing issues with rational points. 

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->


Problem statement:
max bending moment bug
```python
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols

L1 = 6.2
L2 = 1.3

E, I = symbols('E, I')
R1, R2 = symbols('R1, R2')
trave = Beam(L1+L2, E, I)
trave.apply_load(R1, 0, -1)
trave.apply_load(R2, L1, -1)

trave.apply_load(10,0,0)

trave.bc_deflection.append((0, 0))
trave.bc_deflection.append((L1, 0))
#b.bc_deflection.append((L1+L2, 0))

trave.solve_for_reaction_loads(R1, R2)
trave.reaction_loads

trave.plot_shear_force() 
trave.plot_bending_moment() 

# Print max shear force and its location on beam
print(trave.max_shear_force())

# Print max bending moment and its location on beam.
print(trave.max_bmoment())

with this code, max_bmoment gives right maximum moment
if I declare end of distributed load, i can't print maximum moment

from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols

L1 = 6.2
L2 = 1.3

E, I = symbols('E, I')
R1, R2 = symbols('R1, R2')
trave = Beam(L1+L2, E, I)
trave.apply_load(R1, 0, -1)
trave.apply_load(R2, L1, -1)

trave.apply_load(10,0,0, 7.5)

trave.bc_deflection.append((0, 0))
trave.bc_deflection.append((L1, 0))
#b.bc_deflection.append((L1+L2, 0))

trave.solve_for_reaction_loads(R1, R2)
trave.reaction_loads

trave.plot_shear_force() 
trave.plot_bending_moment() 

# Print max shear force and its location on beam
print(trave.max_shear_force())

# Print max bending moment and its location on beam.
print(trave.max_bmoment())
```

Hints from issue discussion:

Commit message: Apply suggestions from code review

Patch to review:
diff --git a/.mailmap b/.mailmap
index 20b4d5efc52e..f64d79255ec8 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1119,6 +1119,7 @@ Prateek Papriwal <papriwalprateek@gmail.com>
 Praveen Sahu <povinsahu@gmail.com> povinsahu1909 <povinsahu@gmail.com>
 Prayush Dawda <35144226+iamprayush@users.noreply.github.com>
 Prempal Singh <prempal.42@gmail.com>
+Prey Patel <patel.prey@iitgn.ac.in>
 Priit Laes <plaes@plaes.org>
 Prince Gupta <codemastercpp@gmail.com> LAPTOP-AS1M2R8B\codem <codemastercpp@gmail.com>
 Prionti Nasir <pdn3628@rit.edu>
diff --git a/sympy/physics/continuum_mechanics/beam.py b/sympy/physics/continuum_mechanics/beam.py
index fc9c436414fd..a9a9bc20cc9b 100644
--- a/sympy/physics/continuum_mechanics/beam.py
+++ b/sympy/physics/continuum_mechanics/beam.py
@@ -934,8 +934,8 @@ def max_shear_force(self):
             if isinstance(term, Mul):
                 term = term.args[-1]    # SingularityFunction in the term
             singularity.append(term.args[1])
-        singularity.sort()
         singularity = list(set(singularity))
+        singularity.sort()
 
         intervals = []    # List of Intervals with discrete value of shear force
         shear_values = []   # List of values of shear force in each interval
@@ -1018,8 +1018,8 @@ def max_bmoment(self):
             if isinstance(term, Mul):
                 term = term.args[-1]    # SingularityFunction in the term
             singularity.append(term.args[1])
-        singularity.sort()
         singularity = list(set(singularity))
+        singularity.sort()
 
         intervals = []    # List of Intervals with discrete value of bending moment
         moment_values = []   # List of values of bending moment in each interval
@@ -1027,7 +1027,10 @@ def max_bmoment(self):
             if s == 0:
                 continue
             try:
-                moment_slope = Piecewise((float("nan"), x<=singularity[i-1]),(self.shear_force().rewrite(Piecewise), x<s), (float("nan"), True))
+                moment_slope = Piecewise(
+                    (float("nan"), x <= singularity[i - 1]),
+                    (self.shear_force().rewrite(Piecewise), x < s),
+                    (float("nan"), True))
                 points = solve(moment_slope, x)
                 val = []
                 for point in points:
@@ -1037,6 +1040,7 @@ def max_bmoment(self):
                 max_moment = max(val)
                 moment_values.append(max_moment)
                 intervals.append(points[val.index(max_moment)])
+
             # If bending moment in a particular Interval has zero or constant
             # slope, then above block gives NotImplementedError as solve
             # can't represent Interval solutions.
diff --git a/sympy/physics/continuum_mechanics/tests/test_beam.py b/sympy/physics/continuum_mechanics/tests/test_beam.py
index 3ee6b044b448..7891a2a8cf30 100644
--- a/sympy/physics/continuum_mechanics/tests/test_beam.py
+++ b/sympy/physics/continuum_mechanics/tests/test_beam.py
@@ -514,7 +514,9 @@ def test_max_shear_force():
     b.apply_load(R2, l, -1)
     b.apply_load(P, 0, 0, end=l)
     b.solve_for_reaction_loads(R1, R2)
-    assert b.max_shear_force() == (0, l*Abs(P)/2)
+    max_shear = b.max_shear_force()
+    assert max_shear[0] == 0
+    assert simplify(max_shear[1] - (l*Abs(P)/2)) == 0
 
 
 def test_max_bmoment():
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 4f761ee4f68e..84f751122c99 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1417,8 +1417,15 @@ def _solve(f, *symbols, **flags):
         result = set()
         if any(e.is_zero for e, c in f.args):
             f = f.simplify()  # failure imminent w/o help
-        for i, (expr, cond) in enumerate(f.args):
-            if expr.is_zero:
+
+        cond = neg = True
+        for i, (expr, cnd) in enumerate(f.args):
+            # the explicit condition for this expr is the current cond
+            # and none of the previous conditions
+            cond = And(neg, cnd).simplify()
+            neg = And(neg, ~cond)
+
+            if expr.is_zero and cond != False:
                 raise NotImplementedError(filldedent('''
                     An expression is already zero when %s.
                     This means that in this *region* the solution
@@ -1427,10 +1434,7 @@ def _solve(f, *symbols, **flags):
                     interval it might be resolved with simplification
                     of the Piecewise conditions.''' % cond))
             candidates = _vsolve(expr, symbol, **flags)
-            # the explicit condition for this expr is the current cond
-            # and none of the previous conditions
-            args = [~c for _, c in f.args[:i]] + [cond]
-            cond = And(*args)
+
             for candidate in candidates:
                 if candidate in result:
                     # an unconditional value was already there


Code review comments:
diff_hunk: @@ -1417,8 +1417,15 @@ def _solve(f, *symbols, **flags):
         result = set()
         if any(e.is_zero for e, c in f.args):
             f = f.simplify()  # failure imminent w/o help
-        for i, (expr, cond) in enumerate(f.args):
-            if expr.is_zero:
+
+        cond = neg = True
+        for i, (expr, cnd) in enumerate(f.args):
+            # the explicit condition for this expr is the current cond
+            # and none of the previous conditions
+            cond = And(neg, cnd).simplify()
+            neg = And(neg, ~cond)
+
+            if expr.is_zero and cond != False:
path: sympy/solvers/solvers.py
text: @user1:
```suggestion
            cond = And(neg, cnd)
            neg = And(neg, ~cond)

            if expr.is_zero and cond.simplify() != False:
```