Differential equations with rational functions as solution

Question

I have some families of nonlinear, first order differential equations. When I try to use DSolve, I usually get a mess (if anything at all) in terms of InverseFunctions and Solve, but I know that my equations always have a rational function as a solution; I just want to know what it is.

If I put a small bound on the degree, I can write the solution out in terms of unknown coefficients, substitute, expand and collect terms, and solve the resulting system of equations (by hand, even).

But clearly this is unsustainable as the degree goes up and even intractable to do by computer. Is there a better way to solve these kinds of differential equation in Mathematica, perhaps somehow telling it the solution will be of this form? For example

DSolve[(x + x^10) == f'[x] (x + f[x]^10), f[x], x]

doesn't even return anything. But $f(x)=x$ is a solution. Similarly,

DSolve[(x + x^2) == f'[x]*(x^2 + f[x]), f[x], x]

returns

Solve[C[1] == (x BesselJ[1, 2 I Sqrt[x - f[x]]] + 
  I BesselJ[2, 2 I Sqrt[x - f[x]]] Sqrt[x - f[x]])/(
    x BesselY[1, -2 I Sqrt[x - f[x]]] - 
     I BesselY[2, -2 I Sqrt[x - f[x]]] Sqrt[x - f[x]]), f[x]]

but $f(x)=x$ is also a solution. In general I would like to handle higher powers of $f'(x)$ as well. But I'd like to figure out simple cases first. Some more examples:

DSolve[ (1 + x^5) == (f'[x])^10*(1 + f[x]^5), f[x], x]

has $f(x) = x$ as a solution, but DSolve gives some mess like:

f[x] -> InverseFunction[
   1/3 #1 (Hypergeometric2F1[1/5, 9/10, 6/5, -#1^5] + 
       2 (1 + #1^5)^(1/10)) &][
  C[1] - 1/3 x (2 (1 + x^5)^(1/10) + 
      Hypergeometric2F1[1/5, 9/10, 6/5, -x^5])]

which I'm sure is accurate, but not so useful to me. Finally,

DSolve[1 + x^3 == (x^7) (f'[x])^2 + x*f[x], f[x], x]

has solution $f(x) = 1/x$, but DSolve doesn't give any output.

Answer 1

I have no idea how robust the following code is. The idea is to search for singular points and then write f as a polynomial times whatever singularities we found. The degree of the polynomial is bounded by looking at the behaviour at infinity. I've included a last step which checks that the function is a solution, to avoid returning too much junk in case the differential equation did not have a rational solution.

First a tool to solve iteratively a system in which the first equation fixes the first variable, the first two fix the first two etc.

ClearAll[specialsolve];
specialsolve::toomanyvars = "Extra vars `1`.  Rules found: `2`.";
specialsolve[eqs_List, {}, rules_List] := 
 If[And @@ (eqs /. 
     rules), {rules}, {}, {}]; (*Returns a set of solutions.*)
specialsolve[{}, vars_List, rules_List] := 
 Message[specialsolve::toomanyvars, vars, rules];
specialsolve[eqs_List, vars_List, rules_List] :=
  Join @@ 
   Map[specialsolve[Rest[eqs], Rest[vars], Join[rules, #]] &, 
    Solve[eqs[[1]] /. rules, vars[[1]]]];

Next, the horrible code.

ClearAll[rationalDSolve];
rationalDSolve[eq_] := 
 Module[{terms, coefs, powers, specialas, valuesofbforgenerica, 
   valuesofbforspecialas, valuesofbanda, valuesint, 
   inversedenominator, powersnearinfinity, maxpower, coefseqs, 
   solutions},
  (*Use a toy Ansatz of f[x]=c(x-a)^
  b to find behaviour near some point x\[Equal]a.*)
  terms = 
   Expand[eq /. f :> (c (# - a)^b &) /. x :> y + a] /. 
    HoldPattern[Plus[t___] | t___] :> {t};
  coefs = Replace[terms, {p_. y^(q_: 1) :> {q, p}, p_ :> {0, p}}, {1}];
  powers = List @@ Simplify[Min[coefs[[All, 1]]]];
  (*We've found the list of posisble leading powers of y=(x-
  a) for generic a (which one in this list is really the leading \
power depends on b.  
  Determine values of a for which the coefficient of one of these \
powers vanishes.*)
  specialas = 
   DeleteDuplicates[
    Flatten[Map[Solve[Coefficient[Plus @@ terms, y, #] == 0, a] &, 
      powers]]];
  (*For generic a, 
  values of b for which two terms turn out to have the same power are \
candidates for cancellations.  TODO: 
  check that the resulting power is indeed the minimum, 
  otherwise these solutions are useless (but not harmful).*)
  valuesofbforgenerica = 
   Flatten[Map[Solve[Equal @@ #, b] &, Subsets[powers, {2}]], {1, 2}];
  (*For special a, same story.  
  There is probably a more efficient way of doing all that.*)
  valuesofbforspecialas = Flatten[
    specialas /. (a -> val_) :>
      With[{terms = 
         Expand[eq /. f :> (c (# - val)^b &) /. x :> y + val] /. 
          HoldPattern[Plus[t___] | t___] :> {t}},
       With[{coefs = 
          Replace[terms, {p_. y^(q_: 1) :> {q, p}, 
            p_ :> {0, p}}, {1}]},
        With[{powers = List @@ Simplify[Min[coefs[[All, 1]]]]},
         Flatten[
          Map[Solve[{a == val, Equal @@ #}, {a, b}] &, 
           Subsets[powers, {2}]], {1, 2}]]]], {1, 2}];
  (*List of lists of rules {...,{a\[Rule]-1,b\[Rule]3/2},{a\[Rule]-I,
  b\[Rule]3/2},...} for possible leading behaviours.*)
  valuesofbanda = Join[valuesofbforgenerica, valuesofbforspecialas];
  (*Let's focus on solutions that are rational functions, 
  so only integer b.  Also we only care about singularities so b<0.*)


  valuesint = 
   Cases[valuesofbanda, {___, b -> n_Integer?Negative, ___}];
  inversedenominator = Product[(x - a)^b /. vals, {vals, valuesint}];
  (*Powers near x\[Rule]Infinity.  
  Now the dominant power is the Max.  
  Directly solve for b integer (we assume rational functions).  
  We care about the maximum possible exponent so take Max afterwards \
too.*)
  powersnearinfinity = List @@ Simplify[Max[coefs[[All, 1]]]];
  maxpower = 
   Max[b /. 
      Flatten[Map[Solve[Equal @@ #, b, Integers] &, 
        Subsets[powersnearinfinity, {2}]], {1, 2}]] + 
    Sum[-b /. val, {val, valuesint}];
  coefseqs = 
   Map[# == 0 &, 
    CoefficientList[
     eq /. f :> (Sum[
            a[i] #^i, {i, 0, maxpower}] inversedenominator /. 
          x -> # &), x]];
  solutions = 
   inversedenominator Sum[a[i] x^i, {i, 0, maxpower}] /. 
    specialsolve[coefseqs, Table[a[i], {i, 0, maxpower}], {}];
  Cases[solutions, 
   sol_ /; FullSimplify[0 == eq /. f :> (sol /. x -> # &)]]]

Let's test.

In[86]:= rationalDSolve[961 x^96 - 899 x^64 - 61 x^32 - 1 - x^5 f'[x]^2 (1 + f[x])]

Out[86]= {(-1 - x + x^32)/x,
  ((-1)^(1/3) - x - (-1)^(1/3) x^32)/x,
  (-(-1)^(2/3) - x + (-1)^(2/3) x^32)/x}

This takes a couple seconds on my desktop. Be patient.