(ode, sol, func=None, order='auto', solve_for_func=True)
r"""
    Substitutes ``sol`` into ``ode`` and checks that the result is ``0``.

    This only works when ``func`` is one function, like `f(x)`.  ``sol`` can
    be a single solution or a list of solutions.  Each solution may be an
    :py:class:`~sympy.core.relational.Equality` that the solution satisfies,
    e.g. ``Eq(f(x), C1), Eq(f(x) + C1, 0)``; or simply an
    :py:class:`~sympy.core.expr.Expr`, e.g. ``f(x) - C1``. In most cases it
    will not be necessary to explicitly identify the function, but if the
    function cannot be inferred from the original equation it can be supplied
    through the ``func`` argument.

    If a sequence of solutions is passed, the same sort of container will be
    used to return the result for each solution.

    It tries the following methods, in order, until it finds zero equivalence:

    1. Substitute the solution for `f` in the original equation.  This only
       works if ``ode`` is solved for `f`.  It will attempt to solve it first
       unless ``solve_for_func == False``.
    2. Take `n` derivatives of the solution, where `n` is the order of
       ``ode``, and check to see if that is equal to the solution.  This only
       works on exact ODEs.
    3. Take the 1st, 2nd, ..., `n`\th derivatives of the solution, each time
       solving for the derivative of `f` of that order (this will always be
       possible because `f` is a linear operator). Then back substitute each
       derivative into ``ode`` in reverse order.

    This function returns a tuple.  The first item in the tuple is ``True`` if
    the substitution results in ``0``, and ``False`` otherwise. The second
    item in the tuple is what the substitution results in.  It should always
    be ``0`` if the first item is ``True``. Note that sometimes this function
    will ``False``, but with an expression that is identically equal to ``0``,
    instead of returning ``True``.  This is because
    :py:meth:`~sympy.simplify.simplify.simplify` cannot reduce the expression
    to ``0``.  If an expression returned by this function vanishes
    identically, then ``sol`` really is a solution to ``ode``.

    If this function seems to hang, it is probably because of a hard
    simplification.

    To use this function to test, test the first item of the tuple.

    Examples
    ========

    >>> from sympy import Eq, Function, checkodesol, symbols
    >>> x, C1 = symbols('x,C1')
    >>> f = Function('f')
    >>> checkodesol(f(x).diff(x), Eq(f(x), C1))
    (True, 0)
    >>> assert checkodesol(f(x).diff(x), C1)[0]
    >>> assert not checkodesol(f(x).diff(x), x)[0]
    >>> checkodesol(f(x).diff(x, 2), x**2)
    (False, 2)
r"""
    Substitutes ``sol`` into ``ode`` and checks that the result is ``0``.
r Substitutes sol into ode and checks that the result is 0 .
def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
    r"""
    Substitutes ``sol`` into ``ode`` and checks that the result is ``0``.

    This only works when ``func`` is one function, like `f(x)`.  ``sol`` can
    be a single solution or a list of solutions.  Each solution may be an
    :py:class:`~sympy.core.relational.Equality` that the solution satisfies,
    e.g. ``Eq(f(x), C1), Eq(f(x) + C1, 0)``; or simply an
    :py:class:`~sympy.core.expr.Expr`, e.g. ``f(x) - C1``. In most cases it
    will not be necessary to explicitly identify the function, but if the
    function cannot be inferred from the original equation it can be supplied
    through the ``func`` argument.

    If a sequence of solutions is passed, the same sort of container will be
    used to return the result for each solution.

    It tries the following methods, in order, until it finds zero equivalence:

    1. Substitute the solution for `f` in the original equation.  This only
       works if ``ode`` is solved for `f`.  It will attempt to solve it first
       unless ``solve_for_func == False``.
    2. Take `n` derivatives of the solution, where `n` is the order of
       ``ode``, and check to see if that is equal to the solution.  This only
       works on exact ODEs.
    3. Take the 1st, 2nd, ..., `n`\th derivatives of the solution, each time
       solving for the derivative of `f` of that order (this will always be
       possible because `f` is a linear operator). Then back substitute each
       derivative into ``ode`` in reverse order.

    This function returns a tuple.  The first item in the tuple is ``True`` if
    the substitution results in ``0``, and ``False`` otherwise. The second
    item in the tuple is what the substitution results in.  It should always
    be ``0`` if the first item is ``True``. Note that sometimes this function
    will ``False``, but with an expression that is identically equal to ``0``,
    instead of returning ``True``.  This is because
    :py:meth:`~sympy.simplify.simplify.simplify` cannot reduce the expression
    to ``0``.  If an expression returned by this function vanishes
    identically, then ``sol`` really is a solution to ``ode``.

    If this function seems to hang, it is probably because of a hard
    simplification.

    To use this function to test, test the first item of the tuple.

    Examples
    ========

    >>> from sympy import Eq, Function, checkodesol, symbols
    >>> x, C1 = symbols('x,C1')
    >>> f = Function('f')
    >>> checkodesol(f(x).diff(x), Eq(f(x), C1))
    (True, 0)
    >>> assert checkodesol(f(x).diff(x), C1)[0]
    >>> assert not checkodesol(f(x).diff(x), x)[0]
    >>> checkodesol(f(x).diff(x, 2), x**2)
    (False, 2)

    """
    if not isinstance(ode, Equality):
        ode = Eq(ode, 0)
    if func is None:
        try:
            _, func = _preprocess(ode.lhs)
        except ValueError:
            funcs = [s.atoms(AppliedUndef) for s in (
                sol if is_sequence(sol, set) else [sol])]
            funcs = reduce(set.union, funcs, set())
            if len(funcs) != 1:
                raise ValueError(
                    'must pass func arg to checkodesol for this case.')
            func = funcs.pop()
    if not isinstance(func, AppliedUndef) or len(func.args) != 1:
        raise ValueError(
            "func must be a function of one variable, not %s" % func)
    if is_sequence(sol, set):
        return type(sol)([checkodesol(ode, i, order=order, solve_for_func=solve_for_func) for i in sol])

    if not isinstance(sol, Equality):
        sol = Eq(func, sol)
    x = func.args[0]
    s = True
    testnum = 0
    if order == 'auto':
        order = ode_order(ode, func)
    if solve_for_func and not (
            sol.lhs == func and not sol.rhs.has(func)) and not (
            sol.rhs == func and not sol.lhs.has(func)):
        try:
            solved = solve(sol, func)
            if not solved:
                raise NotImplementedError
        except NotImplementedError:
            pass
        else:
            if len(solved) == 1:
                result = checkodesol(ode, Eq(func, solved[0]),
                    order=order, solve_for_func=False)
            else:
                result = checkodesol(ode, [Eq(func, t) for t in solved],
                order=order, solve_for_func=False)

            return result

    while s:
        if testnum == 0:
            # First pass, try substituting a solved solution directly into the
            # ODE. This has the highest chance of succeeding.
            ode_diff = ode.lhs - ode.rhs

            if sol.lhs == func:
                s = sub_func_doit(ode_diff, func, sol.rhs)
            elif sol.rhs == func:
                s = sub_func_doit(ode_diff, func, sol.lhs)
            else:
                testnum += 1
                continue
            ss = simplify(s)
            if ss:
                # with the new numer_denom in power.py, if we do a simple
                # expansion then testnum == 0 verifies all solutions.
                s = ss.expand(force=True)
            else:
                s = 0
            testnum += 1
        elif testnum == 1:
            # Second pass. If we cannot substitute f, try seeing if the nth
            # derivative is equal, this will only work for odes that are exact,
            # by definition.
            s = simplify(
                trigsimp(diff(sol.lhs, x, order) - diff(sol.rhs, x, order)) -
                trigsimp(ode.lhs) + trigsimp(ode.rhs))
            # s2 = simplify(
            #     diff(sol.lhs, x, order) - diff(sol.rhs, x, order) - \
            #     ode.lhs + ode.rhs)
            testnum += 1
        elif testnum == 2:
            # Third pass. Try solving for df/dx and substituting that into the
            # ODE. Thanks to Chris Smith for suggesting this method.  Many of
            # the comments below are his too.
            # The method:
            # - Take each of 1..n derivatives of the solution.
            # - Solve each nth derivative for d^(n)f/dx^(n)
            #   (the differential of that order)
            # - Back substitute into the ODE in decreasing order
            #   (i.e., n, n-1, ...)
            # - Check the result for zero equivalence
            if sol.lhs == func and not sol.rhs.has(func):
                diffsols = {0: sol.rhs}
            elif sol.rhs == func and not sol.lhs.has(func):
                diffsols = {0: sol.lhs}
            else:
                diffsols = {}
            sol = sol.lhs - sol.rhs
            for i in range(1, order + 1):
                # Differentiation is a linear operator, so there should always
                # be 1 solution. Nonetheless, we test just to make sure.
                # We only need to solve once.  After that, we automatically
                # have the solution to the differential in the order we want.
                if i == 1:
                    ds = sol.diff(x)
                    try:
                        sdf = solve(ds, func.diff(x, i))
                        if not sdf:
                            raise NotImplementedError
                    except NotImplementedError:
                        testnum += 1
                        break
                    else:
                        diffsols[i] = sdf[0]
                else:
                    # This is what the solution says df/dx should be.
                    diffsols[i] = diffsols[i - 1].diff(x)

            # Make sure the above didn't fail.
            if testnum > 2:
                continue
            else:
                # Substitute it into ODE to check for self consistency.
                lhs, rhs = ode.lhs, ode.rhs
                for i in range(order, -1, -1):
                    if i == 0 and 0 not in diffsols:
                        # We can only substitute f(x) if the solution was
                        # solved for f(x).
                        break
                    lhs = sub_func_doit(lhs, func.diff(x, i), diffsols[i])
                    rhs = sub_func_doit(rhs, func.diff(x, i), diffsols[i])
                    ode_or_bool = Eq(lhs, rhs)
                    ode_or_bool = simplify(ode_or_bool)

                    if isinstance(ode_or_bool, (bool, BooleanAtom)):
                        if ode_or_bool:
                            lhs = rhs = S.Zero
                    else:
                        lhs = ode_or_bool.lhs
                        rhs = ode_or_bool.rhs
                # No sense in overworking simplify -- just prove that the
                # numerator goes to zero
                num = trigsimp((lhs - rhs).as_numer_denom()[0])
                # since solutions are obtained using force=True we test
                # using the same level of assumptions
                ## replace function with dummy so assumptions will work
                _func = Dummy('func')
                num = num.subs(func, _func)
                ## posify the expression
                num, reps = posify(num)
                s = simplify(num).xreplace(reps).xreplace({_func: func})
                testnum += 1
        else:
            break

    if not s:
        return (True, s)
    elif s is True:  # The code above never was able to change s
        raise NotImplementedError("Unable to test if " + str(sol) +
            " is a solution to " + str(ode) + ".")
    else:
        return (False, s)
def checkodesol ( ode , sol , func = None , order = 'auto' , solve_for_func = True ) : if not isinstance ( ode , Equality ) : ode = Eq ( ode , 0 ) if func is None : try : _ , func = _preprocess ( ode . lhs ) except ValueError : funcs = [ s . atoms ( AppliedUndef ) for s in ( sol if is_sequence ( sol , set ) else [ sol ] ) ] funcs = reduce ( set . union , funcs , set ( ) ) if len ( funcs ) != 1 : raise ValueError ( 'must pass func arg to checkodesol for this case.' ) func = funcs . pop ( ) if not isinstance ( func , AppliedUndef ) or len ( func . args ) != 1 : raise ValueError ( "func must be a function of one variable, not %s" % func ) if is_sequence ( sol , set ) : return type ( sol ) ( [ checkodesol ( ode , i , order = order , solve_for_func = solve_for_func ) for i in sol ] ) if not isinstance ( sol , Equality ) : sol = Eq ( func , sol ) x = func . args [ 0 ] s = True testnum = 0 if order == 'auto' : order = ode_order ( ode , func ) if solve_for_func and not ( sol . lhs == func and not sol . rhs . has ( func ) ) and not ( sol . rhs == func and not sol . lhs . has ( func ) ) : try : solved = solve ( sol , func ) if not solved : raise NotImplementedError except NotImplementedError : pass else : if len ( solved ) == 1 : result = checkodesol ( ode , Eq ( func , solved [ 0 ] ) , order = order , solve_for_func = False ) else : result = checkodesol ( ode , [ Eq ( func , t ) for t in solved ] , order = order , solve_for_func = False ) return result while s : if testnum == 0 : # First pass, try substituting a solved solution directly into the # ODE. This has the highest chance of succeeding. ode_diff = ode . lhs - ode . rhs if sol . lhs == func : s = sub_func_doit ( ode_diff , func , sol . rhs ) elif sol . rhs == func : s = sub_func_doit ( ode_diff , func , sol . lhs ) else : testnum += 1 continue ss = simplify ( s ) if ss : # with the new numer_denom in power.py, if we do a simple # expansion then testnum == 0 verifies all solutions. s = ss . expand ( force = True ) else : s = 0 testnum += 1 elif testnum == 1 : # Second pass. If we cannot substitute f, try seeing if the nth # derivative is equal, this will only work for odes that are exact, # by definition. s = simplify ( trigsimp ( diff ( sol . lhs , x , order ) - diff ( sol . rhs , x , order ) ) - trigsimp ( ode . lhs ) + trigsimp ( ode . rhs ) ) # s2 = simplify( #     diff(sol.lhs, x, order) - diff(sol.rhs, x, order) - \ #     ode.lhs + ode.rhs) testnum += 1 elif testnum == 2 : # Third pass. Try solving for df/dx and substituting that into the # ODE. Thanks to Chris Smith for suggesting this method.  Many of # the comments below are his too. # The method: # - Take each of 1..n derivatives of the solution. # - Solve each nth derivative for d^(n)f/dx^(n) #   (the differential of that order) # - Back substitute into the ODE in decreasing order #   (i.e., n, n-1, ...) # - Check the result for zero equivalence if sol . lhs == func and not sol . rhs . has ( func ) : diffsols = { 0 : sol . rhs } elif sol . rhs == func and not sol . lhs . has ( func ) : diffsols = { 0 : sol . lhs } else : diffsols = { } sol = sol . lhs - sol . rhs for i in range ( 1 , order + 1 ) : # Differentiation is a linear operator, so there should always # be 1 solution. Nonetheless, we test just to make sure. # We only need to solve once.  After that, we automatically # have the solution to the differential in the order we want. if i == 1 : ds = sol . diff ( x ) try : sdf = solve ( ds , func . diff ( x , i ) ) if not sdf : raise NotImplementedError except NotImplementedError : testnum += 1 break else : diffsols [ i ] = sdf [ 0 ] else : # This is what the solution says df/dx should be. diffsols [ i ] = diffsols [ i - 1 ] . diff ( x ) # Make sure the above didn't fail. if testnum > 2 : continue else : # Substitute it into ODE to check for self consistency. lhs , rhs = ode . lhs , ode . rhs for i in range ( order , - 1 , - 1 ) : if i == 0 and 0 not in diffsols : # We can only substitute f(x) if the solution was # solved for f(x). break lhs = sub_func_doit ( lhs , func . diff ( x , i ) , diffsols [ i ] ) rhs = sub_func_doit ( rhs , func . diff ( x , i ) , diffsols [ i ] ) ode_or_bool = Eq ( lhs , rhs ) ode_or_bool = simplify ( ode_or_bool ) if isinstance ( ode_or_bool , ( bool , BooleanAtom ) ) : if ode_or_bool : lhs = rhs = S . Zero else : lhs = ode_or_bool . lhs rhs = ode_or_bool . rhs # No sense in overworking simplify -- just prove that the # numerator goes to zero num = trigsimp ( ( lhs - rhs ) . as_numer_denom ( ) [ 0 ] ) # since solutions are obtained using force=True we test # using the same level of assumptions ## replace function with dummy so assumptions will work _func = Dummy ( 'func' ) num = num . subs ( func , _func ) ## posify the expression num , reps = posify ( num ) s = simplify ( num ) . xreplace ( reps ) . xreplace ( { _func : func } ) testnum += 1 else : break if not s : return ( True , s ) elif s is True : # The code above never was able to change s raise NotImplementedError ( "Unable to test if " + str ( sol ) + " is a solution to " + str ( ode ) + "." ) else : return ( False , s )