I'll post a warning first, as what I'm going to do is, for me, experimental: I have just uploaded a maple worksheet (I had to remove the output, otherwise it was apparently too big) and I'm going to try to insert it here. Given my previous failures at these sort of trivial manipulations, I thought I should send a warning that that's what I'm about to try. Fingers crossed...

thanks a lot Robert

yes, the 3-D system can be written:

> f := x-> x^.5-.1*x;
> xd := f(x)-c;
> cd := (diff(f(x),x)-q)*c;
> qd := -(q-.2e-1)*(q-.3e-1)+(q-.25e-1)*(diff(f(x),x)-q)*c/(f(x)-c);

and thus

> qd := -(q-.2e-1)*(q-.3e-1)+(q-.25e-1)*cd/xd;

Therefore, as cd-->0 and xd-->0, we have: cd/xd-->0/0

I must add that I'm interested in non-zero solutions. I am searching for solutions (if any) that satisfy: x>0, c>0, and .2e-1 < q < .3e-1

My conjecture was that the indeterminate form 0/0 could be finite, something like cd/xd-->a, for a some real number to be determined. A solution of {xd=0,cd=0,qd=0} would be a solution of the 4-D system in (x,c,q,a):

c = f(x)
q = f'(x)
a = limit( (f'(x)-q)*c/(f(x)-c) )
a = (q-.2e-1)*(q-.3e-1) / (q-.25e-1)

The limit is to be understood as being taken for (x,c,q) constrained to the xd=0, cd=0, and qd=0 planes.

The plot is not clear enough to settle the matter, but it does look like the 3 planes xd=0, cd=0, and qd=0 do intersect at some points.

I tried a second-order approximation of the indeterminate form (the first-order approx is obviously inconclusive), but what I get is a cubic equation that's quite messy. I would imagine that Maple's numerical solvers would be more powerful than a second-order approximation.

Or are you saying that the above 3-D and/or 4-D systems do not have a (non-zero) solution?

Clearly, if the 0/0 indeterminate form tends to either 0 or infinity, there will be no non-zero solution of the kind I'm looking for (.2e-1 < q < .3e-1). Can 0/0 be finite? Can there be a solution if 0/0 is finite?

many thanks for your help!

I'll be more careful with the html input option. The mess probably happened when I tried to insert html links to the pictures. For some reason that failed. Can you see the picture? located at www. and:

mapleprimes.com/files/9249_pxyz.jpg

I don't seem to have managed to post these plots, I pulled down the insert/link button and pasted the links, but there's obviously more to it than that!

http://www.mapleprimes.com/files/9249_pz1.jpg
http://www.mapleprimes.com/files/9249_pxyz.jpg

When writing up this post, I noticed some odd behavior. My original system is (x,c,q). I thought I should change that to (x,y,z) to be more elegant. And suddenly Maple couldn't solve it anymore...

xx := x^.5-.1*x-y;
yy := (.5/x^.5-.1-z)*y;
zz := -(z-.02)*(z-.03)+(z-.025)*(.5/x^.5-.1-z)*y/(x^.5-.1*x-y);

Digits:=10: interface(displayprecision=10):
ss:= fsolve({xx,yy/y,zz}, {x,y,z}, avoid={x=0, y=0, z=0});
xs:=subs(ss,x): ys:=subs(ss,y): zs:=subs(ss,z):

and now Maple can't find the solution anymore ... or is it me?

Thanks for the explanation jakubi, I must say I never had to rotate figures outside of maple. It would indeed be a great visual help in some cases.

Something not quite right with my attempt to post a link. I'll try this instead:

http://www.mapleprimes.com/files/9249_qdot.ps

Hi jakubi, I'm not sure what you mean, it's probably well over my head -- is it about dynamic plots or something? 

... it depends on what exactly is meant by "manipulation", for simple stuff like adding fancy labels in pstricks, I've exported 3d plots before with postscript, including color, axes, etc. and it looks alright.

Anyway, I'll use this as an opportunity to upload a picture (an experiment):

<a href='http://www.mapleprimes.com/files/9249_qdot.ps'>Download 9249_qdot.ps</a><br/><a href='http://www.mapleprimes.com/viewfile/2966'>View file details</a>

generated with the following options, just in case it can be useful to someone:

axesfonts:=axesfont=[TIMES,ROMAN,6]:
labelfonts:=labelfont=[TIMES,ROMAN,10]:
webblue:=COLOR(RGB,.1,0,.55):
plotopts:=style=patchcontour, shading=none,lightmodel=light3,axes=boxed:

p:=implicitplot3d({f(x,y,z},x=xmin..xmax,y=ymin..ymax,z=zmin..zmax, numpoints=10000, plotopts, colour=webblue): display3d({p}, axesfonts, labelfonts, orientation=[-25,70]);

plotsetup(ps,plotoutput=`qdot.ps`,plotoptions=`color,portrait,noborder,axiswidth=400pt,axisheight=400pt`):
display({p},orientation=[-25,70], projection=0.8,axesfonts,labelfonts,labels=[``,``,``]); plotsetup(default):
 

I have omitted the code for the actual function f(x,y,z), as the above may be adapted to other uses.

P.S. As I preview the message, the inserted picture does not seem to load itself, as it is a postscript file maybe that's normal. I'll post the message anyway.

was looking for that a while back, will be very useful.

Not sure about the formats you are listing here.

I usually export as postscript and manipulate in pstricks or use directly in latex.

something like that:

plotsetup(ps,plotoutput=`myplot_color.ps`,plotoptions=`color,portrait,noborder,axiswidth=400pt,axisheight=400pt`):  display({p},view=[kmin..kmax,cmin..cmax],labels=[``,``],axesfonts,labelfonts); plotsetup(default):

hope this helps, a bit.

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.

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.
 

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

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.

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]);

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]);