Abel ODEs were allegedly solved in 2005, but are not yet available in Maple... any reason?

Dear ODE enthusiasts and Maple developers,

I recently came across the following article which claims to put to rest the quest for solutions to Abel ODEs:

Dimitrios E. Panayotounakos, "Exact analytic solutions of unsolvable classes of first and second order nonlinear ODEs (Part I: Abel’s equations)" http://dx.doi.org/10.1016/j.aml.2004.09.004

I found the method a little difficult to follow, as it involves numerous changes of variables. Has any of you ODE experts ever looked into this method? Does it look like something that would be easily programmed with Maple?

many thanks,

Patrick.

I would like very much

to read this article, but saddly I have not access to it...

Here is a summary of the solution

Here is a summary of the solution described by Panayotounakos in his article:

http://patrick.toche.free.fr/research/papers/abel.pdf

I have not yet written a maple worksheet to check the numerical validity of the reported solution.

Towards a Maple worksheet to check the solution

I'd like to write a worksheet to check the validity (which I trust it must be, as it has been published in a reputable journal) and (above all) ease of implementability with Maple of Panayotounakos's solution of Abel ODEs.

A closed-form solution to the general Abel ODE should be an event celebrated with glee around the world: I'm a little puzzled by the absence of references to Panayotounakos's article, published in 2005 in Applied Mathematics Letters...

Here is a simple worksheet that follows the steps outlined in Panayotounakos, to the best of my understanding. The notation I have adopted below is mostly Panayotounakos's with some small changes. A summary is provided here: http://patrick.toche.free.fr/research/papers/abel.pdf

> restart: with(DEtools): with(plots): plotsetup(default):
> ode:= diff(y(x),x)*y(x)-y(x) = f(x); Abel ODE to be solved following Panoyotounakos's method
> f := x-> x;  Function on rhs of the Abel ODE to be specified. Simple cases are f(x)=0, f(x)=1, f(x)=x, etc
> A := 1:  Constant of integration to be specified. Will depend on (say) the given value of y(0)
> T := x-> ln(x+2*A):     change of variable from x to t
> X := t-> exp(t)-2*A:   reverse-change of variable from t to x
> F := t-> f(X(t)):      writing f as a function of t rather than x
> G := t-> ( ( t*sin(t) + cos(t) ) * Ci(t) + (cos(t))^(2)   ) * ( 4*t*Ci(t) + cos(t) ) / ( 2*(2*t*Ci(t))^(3) ) - 2*F(t); G(t);
> a := -4;
> b := t-> 3 + 4 * (G(t) + F(t)) * exp(-t);
> c := t-> -4 * (G(t) + 2*F(t)) * exp(-t);
> p := t-> b(t) - a^(2)/3;
> q := t-> c(t) - a*b(t)/3 + 2*(a/3)^(3);
> Cubic := r-> r^(3) + p*r + q;
> P := t-> p(t)/3; this normalization makes the cubic-root trigonometric formula cleaner
> Q := t-> q(t)/2; this normalization makes the cubic-root trigonometric formula cleaner
> R:= (P,Q,k)-> 2*(-P)^(1/2) * sin(( arcsin(Q/P^(2)*(-P)^(1/2)) +2*k*pi)/3); trigonometric formula for the cubic roots
> Y := (x,r) -> (1/2) * (x + 2*A) * (r + 1/3);  Panayotounakos formula for the solution y(x) in terms of the cubic-root-function r(x) defined above
> y := (x)-> Y(x,R(P(T(x)),Q(T(x)),0)): y(x); This is/should be the closed-form solution of y(x)
> py:= plot({y(x)}, x=-2..2, color=black): display({py});

 

P.S. I'm the PatrickT who has posted messages here in the past, for some reason my stars vanished after I updated my email...

JacquesC's picture

Contact Will

There aught to be a way to make sure that the account 'PatrickT' keeps its stars properly, even if you change your email address.

On the formula for cubic roots

By the way, I adopted the trigonometric formula for the cubic root, following suggestions found on mapleprimes. The DalFerro-Tartaglia-Cardano formula is given below.

> r:= (P,Q)-> (-Q + (P^(3)+Q^(2))^(1/2) )^(1/3) + (-Q - (P^(3)+Q^(2))^(1/2) )^(1/3); Cubic root in radicals
> R:= (P,Q,K)-> 2*(-P)^(1/2) * sin(( arcsin(Q/P^(2)*(-P)^(1/2)) +2*K*pi)/3); Cubic root in trigonometrics
> p:=-1: q:=-2: Select values of P and Q
> is(p^(3)+q^(2)>0); Check the sign of the Discriminant for given values of P and Q
> k:=0: select the root, k=0,1,2.
> plot1:= plot({r(p,Q)}, Q, color=black): plot2:= plot({R(p,Q,k)}, Q, color=blue):  display({plot1,plot2}, view=[-10..10,-3..3]);

There is a more succinct way

There is a more succinct way of writing the Panayotounakos formula, summarized here:

http://patrick.toche.free.fr/research/papers/abel2.pdf

The following maple worksheet gives the Panayotounakos formula in a succinct way.  I have also plotted the function h(t) that appears as part of the solution.

I still need to compare this with the maple solution of the original ode. It had better be close! Fingers crossed...

> restart: with(DEtools): with(plots): plotsetup(default):
> f := x-> x;  Function on rhs of the Abel ODE to be specified. Simple cases are f(x)=0, f(x)=1, f(x)=x, etc
> A := 1/2:  Constant of integration to be specified. Will depend on (say) the given value of y(0)
> T := x-> ln(x+2*A):   change of variable from x to t
> F := t-> f(exp(t)-2*A): F(t);  Function on rhs of the Abel ODE as a function of t
> R:= (P,Q,k)-> 2*(-P)^(1/2) * sin(( arcsin(Q/P^(2)*(-P)^(1/2)) +2*k*pi)/3); trigonometric formula for the cubic roots
> H := t-> ( ( t*sin(t) + cos(t) ) * Ci(t) + (cos(t))^(2)   ) * ( 4*t*Ci(t) + cos(t) ) / ( 2*(2*t*Ci(t))^(3) );
> P := t-> (1/3)*(-7/3+4*(H(t)-F(t))*exp(-t));
> Q := t-> (1/2)*(-20/27+4/3*(H(t)-4*F(t))*exp(-t));
> Y := (t,r) -> (1/2) * exp(t) * (r + 1/3);
> y := x-> Y(T(x),R(P(T(x)),Q(T(x)),0)): y(x);
> py:= plot({y(x)}, x=0..3, color=blue): display({py});
> pH:= plot({H(t)}, t=0..7, numpoints=1000, color=blue): display({pH}, view=[0.1..6.3,0..10]);

My first approach

As the method is somewhat involved and the notation a bit ambiguous, I wished first to check both the plausibility of this solution and get an understanding of what is involved, before  going to the proof.

Starting with the notations in the Patrick's version, for the equation:

ode:= diff(y(x),x)*y(x)-y(x) = f(x);

the first issue is the change of variables (x,y)->(t,r):

with(PDEtools,dchange):
dchange({x=exp(t)-2*A,y(x)=1/2*exp(t)*(r(t)+1/3)},ode,[t,r],known=f):
collect(4*%,diff,expand@(u->u/exp(t))):
d2:=subs(f(exp(t)-2*A)=f1(t),%);


                      /d      \       2                    4 f1(t)
   d2 := (r(t) + 1/3) |-- r(t)| + r(t)  - 4/3 r(t) - 5/9 = -------
                      \dt     /                            exp(t)

Ie we get another Abel equation for r(t), where f1(t)=f(x(t)). Now, r(t) is given as the real root(s) of the cubic:

eq:=r(t)^3+p(t)*r(t)+q(t);

where p(t) and q(t) are given functions (involving f1(t)), for which I use 'pt' and 'qt':

g:= ( ( t*sin(t) + cos(t) ) * Ci(t) + (cos(t))^(2)   ) * ( 4*t*Ci(t) + cos(t) )*exp(-t) / ( 2*(2*t*Ci(t))^(3) )-2*f1(t);
pt:=3+4*(g+f1(t))*exp(-t)-16/3;
qt:=-4*(g+2*f1(t))*exp(-t)+4/3*(3+4*(g+f1(t))*exp(-t))-2*(4/3)^3;

Generically, there will be one or three families of real solutions (it may depend on t). So, my first check was to choose f1(t)=exp(t), solve for a real root of 'eq' at a "nice" point, use that as IC for 'd2' and compare its numerical solution with the solutions of 'eq' as functions of t.

d2a:=(simplify@subs)(f1(t)=exp(-t),d2);
subs(p(t)=pt,q(t)=qt,f1(t)=exp(-t),eq):
collect(%,r,u->collect(u,exp)):#to make it look nice
eqa:=evalindets(%,`^`,simplify);

= r(t)^3+(-7/3+(1/4*((t*sin(t)+cos(t))*Ci(t)+cos(t)^2)*(4*t*Ci(t)+cos(t))/t^3/Ci(t)^3-4)*exp(-2*t))*r(t)-20/27+(1/12*((t*sin(t)+cos(t))*Ci(t)+cos(t)^2)*(4*t*Ci(t)+cos(t))/t^3/Ci(t)^3-16/3)*exp(-2*t); I choose the IC at t=1

Digits:=14:
eqa1:=eval(eqa,t=1.):
r1q:=fsolve(subs(r(1.)=r1,eqa1),r1);

                        r1q := 1.5164414176292
rn:=dsolve({d2a,r(1.)=r1q},numeric);

And compare the plots of the curves, from 'eqa' with 'implicitplot' and from 'rn' with 'odeplot':

with(plots):
p1:=implicitplot(eqa,t=1..5,r=0..3, gridrefine=2,color=green):
p2:=odeplot(rn,1..5,view=-1..3):
display(p1,p2);

This is not a proof, and there are some discrepancies between the curves. But small discrepancies do not seem to me unexpected due to the numerical errors, mainly at solving the cubic and specially at the switch of its real branch. And that these curves match partially seem to me suggesting that the roots of the cubic are indeed solutions.

discrepancies

Thank you very much jakubi, this is outstanding. Very useful to understand this Panayotounakos formula.  I am also learning a great deal about how to use Maple's capabilities.

The discrepancies between the dsolve,numeric output and the Panayotounakos formula can get quite large, what does it mean?

Take the simple case f(x)=0. In this case, the ODE can be solved in closed form. You'd expect numerical inaccuracies to be minimal in a simple case such as this. Yet, the discrepancies do look rather large.

Consider the case F(t)=0 . I input the Panayotounakos formula a little differently from jakubi, but it amounts to the same: I define a function H(t), in terms of cos(t), sin(t), and Ci(t), and express the Panayotounakos formula in terms of H(t), and F(t). (The formulae could be made even cleaner by defining a function exp(-t)*F(t) and factorizing exp(-2*t), but I haven't bothered) Then I input the trigonometric formula for the cubic roots directly, rather than using implicitplot -- the graph that obtains in this way, for the case F(t)=exp(-t), is the same as the one jakubi has presented. I then looked at different rhs functions, and compared the dsolve,numeric output with the formula -- discrepancies seem to appear even far from the discontinuities.

Starting with the transformed Abel ODE in r(t):

de := (r(t)+1/3)*diff(r(t),t)+r(t)^2-4/3*r(t)-5/9 = 4*F(t)/exp(t); #Transformed Abel ODE

H := t-> ( ( t*sin(t) + cos(t) ) * Ci(t) + (cos(t))^(2)   ) * ( 4*t*Ci(t) + cos(t) ) / ( (t*Ci(t))^(3) ); #Panayotounakos's function
P := (t,H,F) -> -7/9 + 1/12 * H * exp(-2*t) - 4/3 * F * exp(-t) ;
Q := (t,H,F) -> -10/27 + 1/24 * H * exp(-2*t) - 8/3 * F * exp(-t) ;
R := (P,Q,k)-> 2*(-P)^(1/2) * sin(( arcsin(Q/(P)^(2)*(-P)^(1/2)) +2*k*Pi)/3); #trigonometric formula for roots
K := 1: # Selecting one of the real branches of the cubic, K=0,1,2
rsol := (t,F)-> R(P(t,H(t),F),Q(t,H(t),F),K): rsol(t,F): #Panayotounakos's formula for the solution r(t)

The above completes the coding of the Panayotounakos formula for the solution (to the best of my understanding). I now take a function on the rhs of the ODE and compare the formula with dsolve,numeric, following jakubi's suggestions.

FF := t-> 0: #Test function for F(t)

rsol1 := t-> rsol(t,FF(t)): rsol1(t):
r1 := evalf(Re(rsol1(1.)));

                        r1 := 1.16059739086600

de1 := subs( F(t) = FF(t), de );

dsolve(de1); #In some cases the solution can be found explicitly

                r(t) = -1/3, r(t) = 5/3 + exp(-t) _C1

de1sol := dsolve( { de1, r(1.) = r1 }, numeric );

p1 := plot( rsol1(t) , t=0.2..4*Pi, resolution=1000, color=green ):

p2 := odeplot( de1sol, 0.2..4*Pi, resolution=1000, color=blue ):
 

display({p1,p2}, view=[1..4*Pi,0..3]);

 

 

 

Anyhow, I need to think about all this some more. I managed to post a picture, so that's not bad going for one day.

Thanks jakubi for your great help!

Patrick.

what does it mean?

A first guess about it. In my choice f1(t)=exp(-t) [there was a typo above], the rhs tends fast to 0, so it approches to your choice f1(t)=0. Then, for f1(t)=0, the Abel eq. for r(t)

d2 := (r(t)+1/3)*diff(r(t),t)+r(t)^2-4/3*r(t)-5/9 = 4*f1(t)/exp(t);

admits two constant solutions, the roots of a quadratic:

eval(d2,[diff=0,f1=0]);
solve(%,r(t));

                          2
                      r(t)  - 4/3 r(t) - 5/9 = 0

                              5/3, -1/3

So, the issue is that the solution for the Abel eq. for r(t) and a root of the cubic for r(t) are close to 5/3 for these examples when t is large enough. But there is a disagreement where r'(t)<>0 (ie near the initial point). So, something is wrong with the solution in the cubic. May be a typo or an error in the calculations. I have not checked them yet.

PS. Asymptotically, for t-> infinity, this root 5/3 of the quadratic is a root of the cubic. Using the definitions of above:

eqs:=subs(r(t)=5/3,p(t)=pt,q(t)=qt,f1(t)=0,eq);
                                                2
  eqs := 1/2 ((t sin(t) + cos(t)) Ci(t) + cos(t) )

                                    2   /   3      3
        (4 t Ci(t) + cos(t)) exp(-t)   /  (t  Ci(t) )
                                      /

series(eqs,t=infinity,2);
                                O(1)
                               -------
                                     2
                               exp(t)

Thanks JacquesC, I only just

Thanks JacquesC, I only just noticed your message: I have contacted Will actually. I should get the stars back soon. It was only a handful of stars...

Re: Abel ODEs were ellegedly solved in 2005, ...

Hi,
I do not have access to the article you mention, but read your summary, the one you posted in your website as http://patrick.toche.free.fr/research/papers/abel2.pdf

Maple has what I believe are the most advanced algorithms for Abel ODE classes of equations. By most advanced I mean both with respect to other computer algebra systems and as found in the literature. By classes I mean not just one Abel ODE but the family generated from it by changing variables using x-> F(x), y(x) -> ((P+Q*y)/(R+S*y))(x);

If compared with a single ODE, the family involves five additional functions and the algorithms can solve the problem for arbitrary values of them. In practice, when these five functions involve mathematical functions, the algorithm be computationally  too expensive and the gcd and zero recognitions computations with math functions are prone to fail.

Bearing that in mind, Abel equation research splits into two different problems:

1) to find an Abel equation which does not belong to a known solvable class - let's call this ODE "a new seed"

2) to develop an algorithmic way to solve the whole Abel class defined around this seed (ie: to determine the values of the functions F, P, Q, R, S so that a given ODE is reduced to the "new seed"

If you have 1), you can claim that have discovered a new solvable seed. If in addition you have 2), you can fairly claim that a new Abel class is now solvable (in both cases you cannot claim that "Abel ODEs were solved").

Now regarding the ODE example you post,in your pdf, y' y - y = f(x), you do not show anything for 2) and also do not show anything for 1). Regarding 1), it is necessary, first of all, to show a verifiable solution (you have not) - that is a solution in explicit or explicit form that either cancels the equation or that can be used to reconstruct the equation using different approaches; then you need to show that this ODE does not belong to one of the already known solvable classes AIL, AIR, or to one of the solvable AIA subclasses, all of them defined in "Cheb-Terrab, E.S., and Roche, A.D. "An Abel ordinary differential equation class generalizing known integrable classes." European Journal of Applied Mathematics. Vol. 14. (2003): 217-229.".  Years have passed since that paper, that actually report results presented in FOCM in the year 2000, and I have not yet found a "new" seed that does not belong to this class, and algorithm (ie: item 2) above) regarding all these solvable classes already exists in Maple.

Note that your ODE is already AIL or AIR for various values of f(x), for instance f(x) = x, 1/x., 1/x^2 ... (not for x^2), .... So here is the simplification I propose you: can you show how you derive a verifiable solution for f(x) = x^2, that is for y' y - y = x^2, using the method you are presenting? if so we have a good starting point regarding 1).

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Can the formula be made to work?

Hi Edgardo,

Thank you very much for your message. I am aware that Maple has a very advanced algorithm for solving ODEs, and I have benefited from it for some time. I am also aware of your research in this area, and, in particular, the following papers:

E.S. Cheb-Terrab, A.D. Roche, "An Abel ordinary differential equation class generalizing known integrable classes", European Journal of Applied Mathematics, Vol. 14, No. 2, pp. 217-229 (2003).

E.S. Cheb-Terrab, A.D. Roche, "Abel ODEs: Equivalence and Integrable Classes", Computer Physics Communications, 130 (2000) 204

Precisely for these reasons, I was rather surprised by the claims made in:

D. E. Panayotounakos, "Exact analytic solutions of unsolvable classes of first and second order nonlinear ODEs (Part I: Abel’s equations)". Applied Mathematics Letters. Volume 18, Issue 2, February 2005, Pages 155-162.

that "We provide a mathematical technique leading to the construction of exact analytic solutions of the classes of Abel’s nonlinear ordinary differential equations" and even more surprised to find no references to the method on mapleprimes, some three years after the article was published.

Prof. Panayotounakos claims that the Abel ODE can be solved even in cases that do no fall within the integrable classes you identify in your research. In other words, Panayotounakos claims to have both (1) and (2) *for any seed*. An amazing claim, if correct.

My first reaction was indeed to email you (on 29 August) to ask your opinion about Panayotounakos's formula.

I immediately became interested in putting the formula into a Maple worksheet and, thanks to jakubi, the graph above was produced, which shows that (as you note) the Panayotounakos formula *as I have written it* is not exact.

The question that remains is therefore: does the formula contain some typo that could be corrected or is the formula itself incorrect?

If the Panayotounakos formula is correct, it must surely be of great interest to Maple developers. Unfortunately, I have not been able to follow all the steps followed by Prof. Panayotounakos in his derivations, so I am unable to ascertain the validity of the formula.

The step you propose, namely to "show how you derive a verifiable solution for f(x) = x^2, that is for y' y - y = x^2, using the method you are presenting" is indeed the step that has been attempted here, for the even easier case f(x)=0. But as I remarked, the step fails, so I am left wondering: am I misunderstanding Prof. Panayotounakos's claims? am I misunderstanding the formula he presents? I am, unfortunately, unable to reach a conclusion on this matter at this time...

the search for a formula goes on!

many thanks for your suggestions,

Patrick.
 

I don't think so. I gave a

I don't think so. I gave a look at the actual paper now. In the first place, Abel equations of the second kind, in their most general form, they have a numerator of degree 3 in y in the right hand-side of "y' = ...."  , not of degree 2. See the Maple help page ?odeadvisor[Abel2C]. In the first formula in that paper, however, the author presents a numerator of degree 2 and that is key in all what follows. So to start with, this paper discuss only a particular case of Abel equations.

I say that this particularization is key in what follows because, as the author says, his result assumes that any Abel equation can be rewritten as (y'-1) y = f(x), and that is not possible in the general case unless you solve the equation you are trying to rewrite, itself. Or unless one of the coefficients is particularized in such a way that can be expressed in terms of the other ones.

So, understanding that this paper discusses only the particular cases that can be mapped into (y'-1) y = f(x), it is also the case that the transformation mapping `y'` = (f[2]*y^2+f[1]*y+f[0])/(g[1]*y+g[0]) (that is, without cubic term in the numerator of the rhs) into that normal form requires computing the inverse of the integral of the function f[1], which is not always feasible.

Assuming that achieving that normal form in that particular case is possible, or that a RootOf an uncomputed integral in the coefficients of the ODE is satisfactory for the purpose of representing the normal form, it is still the case that the authors of these papers do not show an ODE being solved with their method (their papers do not contain an example), and that the formula derived for the ODE (y'-1) y = f(x) does not verify even for the simplest forms of f(x).
 
 
> If the Panayotounakos formula is correct, it must surely be of great interest to Maple developers
 
Perhaps this is just my perception but I feel that nowadays many people claim many things, while few however show something new. In this context, a claim, without even one example, is not sufficient. So with due respect I would rephrase your sentence the other way around: show one example solved using those formulas that is not already solved by the software working without human intervention, and that would make those formulas of interest for Maple developers.
 
As said in the other post, here is one example, of the form (y'-1) y = f(x) that these authors claim to solve, an example that is not solved by the software without human intervention, the simplest one that is not of the AIR, AIL classes, both fully solvable, and also not of one the known solvable subclasses of AIA: take f(x) = x^2. Another suggestion: instead of trying to discover a mistake in 8 pages of formulas try writing the authors directly, they will most probably reply you with either the solution or explaining why the method does not apply if that is the case.

Edgardo S. Cheb-Terrab
Physics, Maplesoft
 

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Thanks Edgardo

 

The first thing I did was to email Prof. Panayotounakos. Unfortunately, I have received no reply. It's been more than a month now. He is listed as the corresponding author, and his is the only email available. I thought about contacting the journal to relay my questions, would you recommend it?

There is a brief explanation there http://eqworld.ipmnet.ru/en/solutions/ode/ode0124.pdf and there  http://eqworld.ipmnet.ru/en/solutions/ode/ode0125.pdf on how to transform Abel equations of the second kind into the canonical form (y'-1)y = f(x).

Thanks for your suggestions Edgardo.

Patrick.

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}