def _eval_interval(self, x, a, b):
        """
        Returns evaluation over an interval.  For most functions this is:

        self.subs(x, b) - self.subs(x, a),

        possibly using limit() if NaN is returned from subs.

        If b or a is None, it only evaluates -self.subs(x, a) or self.subs(b, x),
        respectively.

        """
        from sympy.series import limit
        if a is None:
            A = 0
        else:
            A = self.subs(x, a)

        if A is S.NaN:
            A = limit(self, x, a)
            if A is S.NaN:
                return self
        if b is None:
            B = 0
        else:
            B = self.subs(x, b)

        if B is S.NaN:
            B = limit(self, x, b)
        if B is S.NaN:
            return self
        return B - A

    def doit(self, **hints):
        if not hints.get('integrals', True):
            return self

        function = self.function

        for x,ab in self.limits:
            antideriv = self._eval_integral(function, x)

            if antideriv is None:
                return self
            else:
                if ab is None:
                    function = antideriv
                else:
                    a,b = ab
                    A = antideriv.subs(x, a)

                    if isinstance(A, Basic.NaN):
                        A = limit(antideriv, x, a)
                    if isinstance(A, Basic.NaN):
                        return self

                    B = antideriv.subs(x, b)

                    if isinstance(B, Basic.NaN):
                        B = limit(antideriv, x, b)
                    if isinstance(B, Basic.NaN):
                        return self

                    function = B - A

        return function

    def max_bmoment(self):
        """Returns maximum Shear force and its coordinate
        in the Beam object."""
        from sympy import solve, Mul, Interval
        bending_curve = self.bending_moment()
        x = self.variable

        terms = bending_curve.args
        singularity = []        # Points at which bending moment changes
        for term in terms:
            if isinstance(term, Mul):
                term = term.args[-1]    # SingularityFunction in the term
            singularity.append(term.args[1])
        singularity.sort()
        singularity = list(set(singularity))

        intervals = []    # List of Intervals with discrete value of bending moment
        moment_values = []   # List of values of bending moment in each interval
        for i, s in enumerate(singularity):
            if s == 0:
                continue
            try:
                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:
                    val.append(bending_curve.subs(x, point))
                points.extend([singularity[i-1], s])
                val.extend([limit(bending_curve, x, singularity[i-1], '+'), limit(bending_curve, x, s, '-')])
                val = list(map(abs, val))
                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.
            except NotImplementedError:
                initial_moment = limit(bending_curve, x, singularity[i-1], '+')
                final_moment = limit(bending_curve, x, s, '-')
                # If bending_curve has a constant slope(it is a line).
                if bending_curve.subs(x, (singularity[i-1] + s)/2) == (initial_moment + final_moment)/2 and initial_moment != final_moment:
                    moment_values.extend([initial_moment, final_moment])
                    intervals.extend([singularity[i-1], s])
                else:    # bending_curve has same value in whole Interval
                    moment_values.append(final_moment)
                    intervals.append(Interval(singularity[i-1], s))

        moment_values = list(map(abs, moment_values))
        maximum_moment = max(moment_values)
        point = intervals[moment_values.index(maximum_moment)]
        return (point, maximum_moment)

 def calc_limit(a, b):
     """replace x with a, using subs if possible, otherwise limit
     where sign of b is considered"""
     wok = inverse_mapping.subs(x, a)
     if not wok is S.NaN:
         return wok
     return limit(sign(b)*inverse_mapping, x, a)

 def calc_limit(a, b):
     """replace x with a, using subs if possible, otherwise limit
     where sign of b is considered"""
     wok = inverse_mapping.subs(x, a)
     if wok is S.NaN or wok.is_bounded is False and a.is_bounded:
         return limit(sign(b)*inverse_mapping, x, a)
     return wok

    def solve_for_reaction_loads(self, *reactions):
        """
        Solves for the reaction forces.

        Examples
        ========
        There is a beam of length 30 meters. A moment of magnitude 120 Nm is
        applied in the clockwise direction at the end of the beam. A pointload
        of magnitude 8 N is applied from the top of the beam at the starting
        point. There are two simple supports below the beam. One at the end
        and another one at a distance of 10 meters from the start. The
        deflection is restricted at both the supports.

        Using the sign convention of upward forces and clockwise moment
        being positive.

        >>> from sympy.physics.continuum_mechanics.beam import Beam
        >>> from sympy import symbols, linsolve, limit
        >>> E, I = symbols('E, I')
        >>> R1, R2 = symbols('R1, R2')
        >>> b = Beam(30, E, I)
        >>> b.apply_load(-8, 0, -1)
        >>> b.apply_load(R1, 10, -1)  # Reaction force at x = 10
        >>> b.apply_load(R2, 30, -1)  # Reaction force at x = 30
        >>> b.apply_load(120, 30, -2)
        >>> b.bc_deflection = [(10, 0), (30, 0)]
        >>> b.load
        R1*SingularityFunction(x, 10, -1) + R2*SingularityFunction(x, 30, -1)
            - 8*SingularityFunction(x, 0, -1) + 120*SingularityFunction(x, 30, -2)
        >>> b.solve_for_reaction_loads(R1, R2)
        >>> b.reaction_loads
        {R1: 6, R2: 2}
        >>> b.load
        -8*SingularityFunction(x, 0, -1) + 6*SingularityFunction(x, 10, -1)
            + 120*SingularityFunction(x, 30, -2) + 2*SingularityFunction(x, 30, -1)
        """
        if self._composite_type == "hinge":
            return self._solve_hinge_beams(*reactions)

        x = self.variable
        l = self.length
        shear_curve = limit(self.shear_force(), x, l)
        moment_curve = limit(self.bending_moment(), x, l)
        reaction_values = linsolve([shear_curve, moment_curve], reactions).args
        self._reaction_loads = dict(zip(reactions, reaction_values[0]))
        self._load = self._load.subs(self._reaction_loads)

 def _calc_limit_1(F, a, b):
     """
     replace d with a, using subs if possible, otherwise limit
     where sign of b is considered
     """
     wok = F.subs(d, a)
     if wok is S.NaN or wok.is_bounded is False and a.is_bounded:
         return limit(sign(b)*F, d, a)
     return wok

    def _eval_imageset(self, f):
        from sympy.functions.elementary.miscellaneous import Min, Max
        from sympy.solvers import solve
        from sympy.core.function import diff
        from sympy.series import limit
        from sympy.calculus.singularities import singularities
        # TODO: handle piecewise defined functions
        # TODO: handle functions with infinitely many solutions (eg, sin, tan)
        # TODO: handle multivariate functions

        expr = f.expr
        if len(expr.free_symbols) > 1 or len(f.variables) != 1:
            return
        var = f.variables[0]

        if not self.start.is_comparable or not self.end.is_comparable:
            return

        try:
            sing = [x for x in singularities(expr, var) if x.is_real and x in self]
        except NotImplementedError:
            return

        if self.left_open:
            _start = limit(expr, var, self.start, dir="+")
        elif self.start not in sing:
            _start = f(self.start)
        if self.right_open:
            _end = limit(expr, var, self.end, dir="-")
        elif self.end not in sing:
            _end = f(self.end)

        if len(sing) == 0:
            solns = solve(diff(expr, var), var)

            extr = [_start, _end] + [f(x) for x in solns
                                     if x.is_real and x in self]
            start, end = Min(*extr), Max(*extr)

            left_open, right_open = False, False
            if _start <= _end:
                # the minimum or maximum value can occur simultaneously
                # on both the edge of the interval and in some interior
                # point
                if start == _start and start not in solns:
                    left_open = self.left_open
                if end == _end and end not in solns:
                    right_open = self.right_open
            else:
                if start == _end and start not in solns:
                    left_open = self.right_open
                if end == _start and end not in solns:
                    right_open = self.left_open

            return Interval(start, end, left_open, right_open)
        else:
            return imageset(f, Interval(self.start, sing[0],
                                        self.left_open, True)) + \
                Union(*[imageset(f, Interval(sing[i], sing[i + 1]), True, True)
                        for i in range(1, len(sing) - 1)]) + \
                imageset(f, Interval(sing[-1], self.end, True, self.right_open))

    def solve_for_reaction_loads(self, *reactions):
        """
        Solves for the reaction forces.

        Examples
        ========
        There is a beam of length 30 meters. A moment of magnitude 120 Nm is
        applied in the clockwise direction at the end of the beam. A pointload
        of magnitude 8 N is applied from the top of the beam at the starting
        point. There are two simple supports below the beam. One at the end
        and another one at a distance of 10 meters from the start. The
        deflection is restricted at both the supports.

        Using the sign convention of upward forces and clockwise moment
        being positive.

        >>> from sympy.physics.continuum_mechanics.beam import Beam
        >>> from sympy import symbols, linsolve, limit
        >>> E, I = symbols('E, I')
        >>> R1, R2 = symbols('R1, R2')
        >>> b = Beam(30, E, I)
        >>> b.apply_load(-8, 0, -1)
        >>> b.apply_load(R1, 10, -1)  # Reaction force at x = 10
        >>> b.apply_load(R2, 30, -1)  # Reaction force at x = 30
        >>> b.apply_load(120, 30, -2)
        >>> b.bc_deflection = [(10, 0), (30, 0)]
        >>> b.load
        R1*SingularityFunction(x, 10, -1) + R2*SingularityFunction(x, 30, -1)
            - 8*SingularityFunction(x, 0, -1) + 120*SingularityFunction(x, 30, -2)
        >>> b.solve_for_reaction_loads(R1, R2)
        >>> b.reaction_loads
        {R1: 6, R2: 2}
        >>> b.load
        -8*SingularityFunction(x, 0, -1) + 6*SingularityFunction(x, 10, -1)
            + 120*SingularityFunction(x, 30, -2) + 2*SingularityFunction(x, 30, -1)
        """
        if self._composite_type == "hinge":
            return self._solve_hinge_beams(*reactions)

        x = self.variable
        l = self.length
        C3 = Symbol('C3')
        C4 = Symbol('C4')

        shear_curve = limit(self.shear_force(), x, l)
        moment_curve = limit(self.bending_moment(), x, l)

        slope_eqs = []
        deflection_eqs = []

        slope_curve = integrate(self.bending_moment(), x) + C3
        for position, value in self._boundary_conditions['slope']:
            eqs = slope_curve.subs(x, position) - value
            slope_eqs.append(eqs)

        deflection_curve = integrate(slope_curve, x) + C4
        for position, value in self._boundary_conditions['deflection']:
            eqs = deflection_curve.subs(x, position) - value
            deflection_eqs.append(eqs)

        solution = list((linsolve([shear_curve, moment_curve] + slope_eqs
                            + deflection_eqs, (C3, C4) + reactions).args)[0])
        solution = solution[2:]

        self._reaction_loads = dict(zip(reactions, solution))
        self._load = self._load.subs(self._reaction_loads)

    def _solve_hinge_beams(self, *reactions):
        """Method to find integration constants and reactional variables in a
        composite beam connected via hinge.
        This method resolves the composite Beam into its sub-beams and then
        equations of shear force, bending moment, slope and deflection are
        evaluated for both of them separately. These equations are then solved
        for unknown reactions and integration constants using the boundary
        conditions applied on the Beam. Equal deflection of both sub-beams
        at the hinge joint gives us another equation to solve the system.

        Examples
        ========
        A combined beam, with constant fkexural rigidity E*I, is formed by joining
        a Beam of length 2*l to the right of another Beam of length l. The whole beam
        is fixed at both of its both end. A point load of magnitude P is also applied
        from the top at a distance of 2*l from starting point.

        >>> from sympy.physics.continuum_mechanics.beam import Beam
        >>> from sympy import symbols
        >>> E, I = symbols('E, I')
        >>> l=symbols('l', positive=True)
        >>> b1=Beam(l ,E,I)
        >>> b2=Beam(2*l ,E,I)
        >>> b=b1.join(b2,"hinge")
        >>> M1, A1, M2, A2, P = symbols('M1 A1 M2 A2 P')
        >>> b.apply_load(A1,0,-1)
        >>> b.apply_load(M1,0,-2)
        >>> b.apply_load(P,2*l,-1)
        >>> b.apply_load(A2,3*l,-1)
        >>> b.apply_load(M2,3*l,-2)
        >>> b.bc_slope=[(0,0), (3*l, 0)]
        >>> b.bc_deflection=[(0,0), (3*l, 0)]
        >>> b.solve_for_reaction_loads(M1, A1, M2, A2)
        >>> b.reaction_loads
        {A1: -5*P/18, A2: -13*P/18, M1: 5*P*l/18, M2: -4*P*l/9}
        >>> b.slope()
        Piecewise(((5*P*l*SingularityFunction(x, 0, 1)/18 - 5*P*SingularityFunction(x, 0, 2)/36
        + 5*P*SingularityFunction(x, l, 2)/36)/(E*I), l >= x), ((P*l**2/18 - 4*P*l*SingularityFunction(-l +
        x, 2*l, 1)/9 - 5*P*SingularityFunction(-l + x, 0, 2)/36 + P*SingularityFunction(-l + x, l, 2)/2
        - 13*P*SingularityFunction(-l + x, 2*l, 2)/36)/(E*I), x < 3*l))
        >>> b.deflection()
        Piecewise(((5*P*l*SingularityFunction(x, 0, 2)/36 - 5*P*SingularityFunction(x, 0, 3)/108
        + 5*P*SingularityFunction(x, l, 3)/108)/(E*I), l >= x), ((5*P*l**3/54 + P*l**2*(-l + x)/18
        - 2*P*l*SingularityFunction(-l + x, 2*l, 2)/9 - 5*P*SingularityFunction(-l + x, 0, 3)/108
        + P*SingularityFunction(-l + x, l, 3)/6 - 13*P*SingularityFunction(-l + x, 2*l, 3)/108)/(E*I), x < 3*l))
        """
        x = self.variable
        l = self._hinge_position
        E = self._elastic_modulus
        I = self._second_moment

        if isinstance(I, Piecewise):
            I1 = I.args[0][0]
            I2 = I.args[1][0]
        else:
            I1 = I2 = I

        load_1 = 0       # Load equation on first segment of composite beam
        load_2 = 0       # Load equation on second segment of composite beam

        # Distributing load on both segments
        for load in self.applied_loads:
            if load[1] < l:
                load_1 += load[0]*SingularityFunction(x, load[1], load[2])
                if load[2] == 0:
                    load_1 -= load[0]*SingularityFunction(x, load[3], load[2])
                elif load[2] > 0:
                    load_1 -= load[0]*SingularityFunction(x, load[3], load[2]) + load[0]*SingularityFunction(x, load[3], 0)
            elif load[1] == l:
                load_1 += load[0]*SingularityFunction(x, load[1], load[2])
                load_2 += load[0]*SingularityFunction(x, load[1] - l, load[2])
            elif load[1] > l:
                load_2 += load[0]*SingularityFunction(x, load[1] - l, load[2])
                if load[2] == 0:
                    load_2 -= load[0]*SingularityFunction(x, load[3] - l, load[2])
                elif load[2] > 0:
                    load_2 -= load[0]*SingularityFunction(x, load[3] - l, load[2]) + load[0]*SingularityFunction(x, load[3] - l, 0)

        h = Symbol('h')     # Force due to hinge
        load_1 += h*SingularityFunction(x, l, -1)
        load_2 -= h*SingularityFunction(x, 0, -1)

        eq = []
        shear_1 = integrate(load_1, x)
        shear_curve_1 = limit(shear_1, x, l)
        eq.append(shear_curve_1)
        bending_1 = integrate(shear_1, x)
        moment_curve_1 = limit(bending_1, x, l)
        eq.append(moment_curve_1)

        shear_2 = integrate(load_2, x)
        shear_curve_2 = limit(shear_2, x, self.length - l)
        eq.append(shear_curve_2)
        bending_2 = integrate(shear_2, x)
        moment_curve_2 = limit(bending_2, x, self.length - l)
        eq.append(moment_curve_2)

        C1 = Symbol('C1')
        C2 = Symbol('C2')
        C3 = Symbol('C3')
#.........这里部分代码省略.........

def _set_function(f, x):
    from sympy.functions.elementary.miscellaneous import Min, Max
    from sympy.solvers.solveset import solveset
    from sympy.core.function import diff, Lambda
    from sympy.series import limit
    from sympy.calculus.singularities import singularities
    from sympy.sets import Complement
    # TODO: handle functions with infinitely many solutions (eg, sin, tan)
    # TODO: handle multivariate functions

    expr = f.expr
    if len(expr.free_symbols) > 1 or len(f.variables) != 1:
        return
    var = f.variables[0]

    if expr.is_Piecewise:
        result = S.EmptySet
        domain_set = x
        for (p_expr, p_cond) in expr.args:
            if p_cond is true:
                intrvl = domain_set
            else:
                intrvl = p_cond.as_set()
                intrvl = Intersection(domain_set, intrvl)

            if p_expr.is_Number:
                image = FiniteSet(p_expr)
            else:
                image = imageset(Lambda(var, p_expr), intrvl)
            result = Union(result, image)

            # remove the part which has been `imaged`
            domain_set = Complement(domain_set, intrvl)
            if domain_set.is_EmptySet:
                break
        return result

    if not x.start.is_comparable or not x.end.is_comparable:
        return

    try:
        sing = [i for i in singularities(expr, var)
            if i.is_real and i in x]
    except NotImplementedError:
        return

    if x.left_open:
        _start = limit(expr, var, x.start, dir="+")
    elif x.start not in sing:
        _start = f(x.start)
    if x.right_open:
        _end = limit(expr, var, x.end, dir="-")
    elif x.end not in sing:
        _end = f(x.end)

    if len(sing) == 0:
        solns = list(solveset(diff(expr, var), var))

        extr = [_start, _end] + [f(i) for i in solns
                                 if i.is_real and i in x]
        start, end = Min(*extr), Max(*extr)

        left_open, right_open = False, False
        if _start <= _end:
            # the minimum or maximum value can occur simultaneously
            # on both the edge of the interval and in some interior
            # point
            if start == _start and start not in solns:
                left_open = x.left_open
            if end == _end and end not in solns:
                right_open = x.right_open
        else:
            if start == _end and start not in solns:
                left_open = x.right_open
            if end == _start and end not in solns:
                right_open = x.left_open

        return Interval(start, end, left_open, right_open)
    else:
        return imageset(f, Interval(x.start, sing[0],
                                    x.left_open, True)) + \
            Union(*[imageset(f, Interval(sing[i], sing[i + 1], True, True))
                    for i in range(0, len(sing) - 1)]) + \
            imageset(f, Interval(sing[-1], x.end, True, x.right_open))

 def limit(self, x, xlim, direction='+'):
     """ Compute limit x->xlim.
     """
     from sympy.series.limits import limit
     return limit(self, x, xlim, direction)

    def _eval_imageset(self, f):
        from sympy import Dummy
        from sympy.functions.elementary.miscellaneous import Min, Max
        from sympy.solvers import solve
        from sympy.core.function import diff
        from sympy.series import limit
        from sympy.calculus.singularities import singularities
        # TODO: handle piecewise defined functions
        # TODO: handle functions with infinitely many solutions (eg, sin, tan)
        # TODO: handle multivariate functions

        # var and expr are being defined this way to
        # support Python lambda and not just sympy Lambda
        try:
            var = Dummy()
            expr = f(var)
            if len(expr.free_symbols) > 1:
                raise TypeError
        except TypeError:
            raise NotImplementedError("Sorry, Multivariate imagesets are"
                                      " not yet implemented, you are welcome"
                                      " to add this feature in Sympy")

        if not self.start.is_comparable or not self.end.is_comparable:
            raise NotImplementedError("Sets with non comparable/variable"
                                      " arguments are not supported")

        sing = [x for x in singularities(expr, var) if x.is_real and x in self]

        if self.left_open:
            _start = limit(expr, var, self.start, dir="+")
        elif self.start not in sing:
            _start = f(self.start)
        if self.right_open:
            _end = limit(expr, var, self.end, dir="-")
        elif self.end not in sing:
            _end = f(self.end)

        if len(sing) == 0:
            solns = solve(diff(expr, var), var)

            extr = [_start, _end] + [f(x) for x in solns
                                     if x.is_real and x in self]
            start, end = Min(*extr), Max(*extr)

            left_open, right_open = False, False
            if _start <= _end:
                # the minimum or maximum value can occur simultaneously
                # on both the edge of the interval and in some interior
                # point
                if start == _start and start not in solns:
                    left_open = self.left_open
                if end == _end and end not in solns:
                    right_open = self.right_open
            else:
                if start == _end and start not in solns:
                    left_open = self.right_open
                if end == _start and end not in solns:
                    right_open = self.left_open

            return Interval(start, end, left_open, right_open)
        else:
            return imageset(f, Interval(self.start, sing[0],
                                        self.left_open, True)) + \
                Union(*[imageset(f, Interval(sing[i], sing[i + 1]), True, True)
                        for i in range(1, len(sing) - 1)]) + \
                imageset(f, Interval(sing[-1], self.end, True, self.right_open))