# Implementation of the general solution of Abel differential equation of the second kind in Mathematica

## Introduction

An Abel equation of the second kind in its canonical form is writen as

$$y(x)y'(x)-y(x)=f(x)\quad (6)$$ for arbitrary $$f(x)$$. Despite the fact that this equation has general solutions 1 2, Mathematica's tools, ex. DSolve, do not provide a solution. I tried to implement a code that would solve equations of this type and compare it with the known case of $$f(x)=x$$. That is, I am trying to solve:

$$y(x)y'(x)-y(x)=x$$

following the reference 1 and benchmarking it with the known solution returned with DSOlve. However, I wasn't able to do it. Here is my attempt:

## Code

(*Solving for y*)
ζ[x_] = Log[Abs[x + 2 λ]]
F[x_] = ζ[x]
G[ζ_] =
1/16 ((ζ Sin[ζ] + Cos[ζ]) CosIntegral[ζ] +
Cos[ζ]^2)*(4 ζ CosIntegral[ζ] +
Cos[ζ] )/(ζ CosIntegral[ζ])^3 Exp[-ζ] -
2 F[ζ]

a = -4;
b[ζ_] = 3 + 4 (G[ζ] + F[ζ]) Exp[-ζ]
c[ζ_] = - 4 (G[ζ] + 2 F[ζ]) Exp[-ζ]

p[ζ_] = -a^2/3 + b[ζ]
q[ζ_] = 2 (a/3)^3 - a b[ζ]/3 + c[ζ]
NN = Solve[Z^3 + p Z + q == 0, Z] /. {p -> p[ζ],
q -> q[ζ]} /. ζ -> ζ[x];
y = Table[
Table[1/2 (x + 2 λ) ((Z /. NN[[ii]]) + 1/3), {ii, 1,
3}], {λ, 1, 2}];

(*Benchmarking*)
DSolve[Y'[X] Y[X] - Y[X] == X, Y[X], X]

p1 = Plot[y, {x, 0, 10}]
p2 = Table[
ContourPlot[
1/10 ((5 + Sqrt[5]) Log[
1 + Sqrt[5] - (2 Y)/X] - (-5 + Sqrt[5]) Log[-1 + Sqrt[5] + (
2 Y)/X]) == II - Log[X], {X, 0, 10}, {Y, -10, 15}], {II, 0,
3, 0.2}];
Show[p2]


## Plots

The plots look not alike

## References

1. Dimitrios E. Panayotounakos, "Exact analytic solutions of unsolvable classes of first and second" order nonlinear ODEs (Part I: Abel’s equations)
2. Panayotounakos, Dimitrios E.; Zarmpoutis, Theodoros I, "Construction of Exact Parametric or Closed Form Solutions of Some Unsolvable Classes of Nonlinear ODEs (Abel's Nonlinear ODEs of the First Kind and Relative Degenerate Equations)".
• Mathematica is not a magic box that'll spit out a solution to any problem.Maple cannot do it, either.All computer algebra systems, including Mathematica, are limited in their capabilities. Jul 16 '20 at 15:38
• what is actually the solution you obtained? Hard to read all this code. Can you write down, in mathematics, what is the solution you are comparing to Mathematica? You can post it in Latex at top to make it easier to see. I just solved this ODE. It is easy to solve by hand, since it is homogeneous ode. Jul 16 '20 at 17:57
• @MariuszIwaniuk The question posed by the OP is not whether DSolve can solve the ODE (it can, although implicitly) but whether the solution is correct. Jul 17 '20 at 5:04
• I included the analytical solution valid for arbitrary $f(x)$. bbgodfrey pointed out that DSolve doesn't provide the full solution to the case $f(x)=x$. Now, the problem remains to implement the solution provided in the question. Solving by hand is easy for the case $f(x)=x$ but not in the general case... Jul 17 '20 at 10:59
• As a side note, I find it surprising that Mathematica does not provide solutions to such ODE's Jul 17 '20 at 11:01

DSolve easily provides a solution to the ODE, although the solution is implicit.

(DSolve[y[x] (y'[x] - 1) == x, y, x] /. C[1] -> c)[[1]]
(* 1/10 ((5 + Sqrt[5]) Log[1 + Sqrt[5] - (2 y[x])/x] -
(-5 + Sqrt[5]) Log[-1 + Sqrt[5] + (2 y[x])/x]) == c - Log[x] *)


and that solution, in turn, provides the second plot in the Question. The issue raised by the OP is whether that solution is correct. As shown below, the solution is correct. To begin, it is convenient to simplify (by hand) the solution above to

(* 1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] -
(-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) == c *)


which yields the identical plot. However, solutions also exist in the blank regions of the second plot. Replacing Log[x (1 + Sqrt[5]) - 2 y] by the equivalent I Pi + Log[-x (1 + Sqrt[5]) + 2 y] and absorbing the complex constant into c provides solutions in the upper left of the plot, and instead replacing Log[x (-1 + Sqrt[5]) + 2 y]) by I Pi + Log[-x (-1 + Sqrt[5]) - 2 y] and absorbing the complex constant into c provides solutions in the lower left of the plot. Plotting all these representations of the original solution yields

Show[ContourPlot[Evaluate@Table[1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] -
(-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) == c, {c, 0, 3, 0.5}], {x, 0, 10},
{y, -10, 15}, PlotPoints -> 50, ContourStyle -> Blue],
ContourPlot[Evaluate@Table[1/10 ((5 + Sqrt[5]) Log[-x (1 + Sqrt[5]) + 2 y] -
(-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) == c, {c, 0, 3, .5}], {x, 0, 10},
{y, -10, 15}, PlotPoints -> 50, ContourStyle -> Blue],
ContourPlot[Evaluate@Table[1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] -
(-5 + Sqrt[5]) Log[-x (-1 + Sqrt[5]) - 2 y]) == c, {c, 0, 3, .5}], {x, 0, 10},
{y, -10, 15}, PlotPoints -> 50, ContourStyle -> Blue],
ImageSize -> Large, LabelStyle -> {Black, Bold, 15}]


(Noisiness near x (-1 + Sqrt[5]) + 2 y == 0 may be due to a branch cut there.) Now, the question is whether the solution provided by DSolve is correct. Evidence that it is can be obtained by solving the ODE numerically. For the equivalent of c = 1.5,

y /. FindRoot[(1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] -
(-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) == c) /. {x -> 10, c -> 1.5}, {y, 15}];
Quiet@NDSolveValue[{y[x] (y'[x] - 1) == x, y[10] == %}, y, {x, 0, 10}];
%["Domain"] // Flatten;
pa = Plot[%%[x], Evaluate@Join[{x}, %], PlotRange -> {{0, 10}, {-10, 15}}];

y /. FindRoot[(1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] -
(-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) == c) /. {x -> 10, c -> 1.5}, {y, -6}];
Quiet@NDSolveValue[{y[x] (y'[x] - 1) == x, y[10] == %}, y, {x, 0, 10}];
%["Domain"] // Flatten;
pb = Plot[%%[x], Evaluate@Join[{x}, %], PlotRange -> {{0, 10}, {-10, 15}}];

y /. FindRoot[(1/10 ((5 + Sqrt[5]) Log[-x (1 + Sqrt[5]) + 2 y] -
(-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) == c) /. {x -> 10, c -> 1.5}, {y, 18}];
Quiet@NDSolveValue[{y[x] (y'[x] - 1) == x, y[10] == %}, y, {x, 0, 10}];
%["Domain"] // Flatten;
pc = Plot[%%[x], Evaluate@Join[{x}, %], PlotRange -> {{0, 10}, {-10, 15}}];

y /. FindRoot[(1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] -
(-5 + Sqrt[5]) Log[-x (-1 + Sqrt[5]) - 2 y]) == c) /. {x -> 10, c -> 1.5}, {y, -8}];
Quiet@NDSolveValue[{y[x] (y'[x] - 1) == x, y[10] == %}, y, {x, 0, 10}];
%["Domain"] // Flatten;
pd = Plot[%%[x], Evaluate@Join[{x}, %], PlotRange -> {{0, 10}, {-10, 15}}];

Show[pa, pb, pc, pd, ImageSize -> Large, LabelStyle -> {Black, Bold, 15}]


agreeing well with the c = 1.5 curves of the first plot. I cannot speak to the accuracy of either the analytical solution cited in the question or the numerical implementation of it, except to note that NDSolve does not reproduce the first plot in the Question.

• You can also plot the real part of the constant C[1] (or c) in the contour plot: ContourPlot[ 1/10 ((5 + Sqrt[5]) Log[x (1 + Sqrt[5]) - 2 y] - (-5 + Sqrt[5]) Log[x (-1 + Sqrt[5]) + 2 y]) // Re, {x, -12, 12}, {y, -12, 12}] -- Also the case f[x] = x is a homogeneous equation, so the solvability of this special case should seem pretty obvious to DE folks. Jul 17 '20 at 13:27
• @MichaelE2 Certainly, more compact code. Thanks. Jul 17 '20 at 13:32