# Differential equations with rational functions as solution

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.

• Difficult to say much without a specific Mathematica input. Feb 10, 2016 at 23:10
• Thanks. I added two examples, which are more simple than what I'm trying to solve more generally. But I would be happy if I could even get those solved. Feb 10, 2016 at 23:34
• I'm mostly interested in if there is some way to 'hint' to Mathematica (or some other package) that I only want rational functions as solutions. Feb 11, 2016 at 5:08
• Write the solution out in terms of unknown coefficients, substitute, expand and collect terms, and solve the resulting system of equations is the only approach I can think out, and I think it should work when the degree isn't too high. Can you show an example that is hard to handle with this approach? Feb 11, 2016 at 8:21
• Maybe I'm just not doing it the right way, then. For example, $961 x^{96} - 899 x^{64} - 61 x^{32} -1 = (x^5)*f'(x)^2 * (1 + f(x))$ has solution $(x^{32}-x-1)/x$. But if all I know is that the solution is degree 32, then I have 64 unknowns, and running Expand, Collect, Coefficient, Numerator, Simplify, etc on that expression shows no sign of stopping. It's a lot of terms when it's expanded out. I suppose I could try to let it run overnight. Feb 11, 2016 at 9:27

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,
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.

• Thank you so much, this is very much in the direction of what I was looking for. The fact that it can handle the higher degree cases is very useful. I will have to experiment more with this and see what it can do. Jun 9, 2016 at 22:53