       Re: Help with an ODE

• To: mathgroup at smc.vnet.net
• Subject: [mg111304] Re: Help with an ODE
• From: Leonid Shifrin <lshifr at gmail.com>
• Date: Tue, 27 Jul 2010 04:15:59 -0400 (EDT)

Sam,

My guess is that your equation does not have a solution with your boundary
condition. If you use NDSolve, it will bark on derivative at zero being
infinite.
What I can do is to show you how for example your equation can be solved
for y=1, y'=1 and a=1. We will get an analytical solution for x =
x[y], but
will have to solve it numerically for y = y[x].

There are a number of steps. This is your equation

In:= eq = y''[x]==(1+y'[x]^2+ a y[x]^4)/y[x]

Out= (y^\[Prime]\[Prime])[x]==(1+a y[x]^4+(y^\[Prime])[x]^2)/y[x]

1. First, equation does not depend on x explicitly, so you
can make a substitution y'[x]  = f[y]. Then, y''[x] = 1/2 d/dy f[y]^2.
Making another
substitution f[y]^2 = phi[y], you get:

In:= eq1=eq/.{y'[x]->f[y],y''[x]:> 1/2D[f[y]^2,y],y[x]-> y}

Out= f[y] (f^\[Prime])[y]==(1+a y^4+f[y]^2)/y

2. Setting now g[y] == f^2[y] you will have:

In:= eq2 =Apart/@(eq1/.{f[y]
(f^\[Prime])[y]->1/2D[g[y],y],f[y]^2->g[y]})

Out= (g^\[Prime])[y]/2==(1+a y^4)/y+g[y]/y

3. It is easy to solve, with the solution being:

In:= soly = a y^4 + C y^2 - 1

Out= -1 + a y^4 + y^2 C

where C is some integration constant. You can check that it is indeed a
solution:

In:= eq2 /. {g[y] -> soly, g'[y] -> D[soly, y]} // Simplify

Out= True

4. Now you recall that g[y]  = (dy/dx)^2, and therefore:

In:= eq3 =y'[x]==Sqrt[soly]

Out= (y^\[Prime])[x]==Sqrt[-1+a y^4+y^2 C]

(actually plus or minus Sqrt, and later I will use the minus sign).

5. This is the equation which gives then x as a function of y.

In:= HoldForm[-x+C==-Integrate[Sqrt[a t^4+Ct^2-1],{t,y,Infinity}]]
Out= (* Suppressed *)

The reason that I did it in this way rather than from 0 to y with a plus
sign is
that y[x] will never be zero.

6. The integral can be computed:

In:= int = Integrate[1/Sqrt[a t^4+C t^2-1],t]

Out= -((I Sqrt[1+(2 a t^2)/(C-Sqrt[4 a+C^2])] Sqrt[1+(2 a
t^2)/(C+Sqrt[4 a+C^2])] EllipticF[I ArcSinh[Sqrt t
Sqrt[a/(C+Sqrt[4 a+C^2])]],(C+Sqrt[4 a+C^2])/(C-Sqrt[4
a+C^2])])/(Sqrt Sqrt[-1+a t^4+t^2 C] Sqrt[a/(C+Sqrt[4
a+C^2])]))

7. Integral at infinity can be computed as a series expansion around y = 1/q
when q->0:

In:=
Clear[q];
intinf = Normal[
Assuming[q > 0, Simplify[Series[int /. t -> 1/q, {q, 0, 0}]]]]

Out= (Sqrt Sqrt[a/(
C - Sqrt[
4 a + C^2])] (EllipticK[(
2 Sqrt[4 a + C^2])/(-C + Sqrt[
4 a + C^2])] + ((-I a +
Sqrt[-(a^2/(C + Sqrt[4 a + C^2])^2)] (C + Sqrt[
4 a + C^2])) (-1 +
Sqrt[(C - Sqrt[4 a + C^2])/(C + Sqrt[4 a + C^2])]
Sqrt[(C + Sqrt[4 a + C^2])/(
C - Sqrt[4 a + C^2])]) EllipticK[(
C + Sqrt[4 a + C^2])/(C - Sqrt[4 a + C^2])])/(
2 a)))/Sqrt[a]

8. Integral on the lower limit you get just substituting t->y:

In:= inty = int /. t -> y

Out= -((
I Sqrt[1 + (2 a y^2)/(C - Sqrt[4 a + C^2])] Sqrt[
1 + (2 a y^2)/(C + Sqrt[4 a + C^2])]
EllipticF[
I ArcSinh[Sqrt y Sqrt[a/(C + Sqrt[4 a + C^2])]], (
C + Sqrt[4 a + C^2])/(C - Sqrt[4 a + C^2])])/(
Sqrt Sqrt[-1 + a y^4 + y^2 C] Sqrt[a/(
C + Sqrt[4 a + C^2])]))

9. The final result for the integral then:

In:= intfull  = intinf - inty
Out= (I Sqrt[1+(2 a y^2)/(C-Sqrt[4 a+C^2])] Sqrt[1+(2 a
y^2)/(C+Sqrt[4 a+C^2])] EllipticF[I ArcSinh[Sqrt y
Sqrt[a/(C+Sqrt[4 a+C^2])]],(C+Sqrt[4 a+C^2])/(C-Sqrt[4
a+C^2])])/(Sqrt Sqrt[-1+a y^4+y^2 C] Sqrt[a/(C+Sqrt[4
a+C^2])])+(Sqrt Sqrt[a/(C-Sqrt[4 a+C^2])] (EllipticK[(2 Sqrt[4
a+C^2])/(-C+Sqrt[4 a+C^2])]+((-I a+Sqrt[-(a^2/(C+Sqrt[4
a+C^2])^2)] (C+Sqrt[4 a+C^2])) (-1+Sqrt[(C-Sqrt[4
a+C^2])/(C+Sqrt[4 a+C^2])] Sqrt[(C+Sqrt[4
a+C^2])/(C-Sqrt[4 a+C^2])]) EllipticK[(C+Sqrt[4
a+C^2])/(C-Sqrt[4 a+C^2])])/(2 a)))/Sqrt[a]

10. For a = 1, condition y' = 1 gives (since y = 1) :

In:= 1 == Sqrt[-1 + a + C]

Out= 1 == Sqrt[-1 + a + C]

which gives C -> 1 for a == 1.

11. We have then :

In:= -x+C==intfull

Out= -x+C==(I Sqrt[1+(2 a y^2)/(C-Sqrt[4 a+C^2])] Sqrt[1+(2 a
y^2)/(C+Sqrt[4 a+C^2])] EllipticF[I ArcSinh[Sqrt y
Sqrt[a/(C+Sqrt[4 a+C^2])]],(C+Sqrt[4 a+C^2])/(C-Sqrt[4
a+C^2])])/(Sqrt Sqrt[-1+a y^4+y^2 C] Sqrt[a/(C+Sqrt[4
a+C^2])])+(Sqrt Sqrt[a/(C-Sqrt[4 a+C^2])] (EllipticK[(2 Sqrt[4
a+C^2])/(-C+Sqrt[4 a+C^2])]+((-I a+Sqrt[-(a^2/(C+Sqrt[4
a+C^2])^2)] (C+Sqrt[4 a+C^2])) (-1+Sqrt[(C-Sqrt[4
a+C^2])/(C+Sqrt[4 a+C^2])] Sqrt[(C+Sqrt[4
a+C^2])/(C-Sqrt[4 a+C^2])]) EllipticK[(C+Sqrt[4
a+C^2])/(C-Sqrt[4 a+C^2])])/(2 a)))/Sqrt[a]

12.  At x == 0, and a == 1, y == 1, and we get

In:=
c2rule = C-> intfull/.{y->1,C->1,a->1}

Out=
C->-Sqrt[1/2 (1+Sqrt) (-1-2/(1-Sqrt)) (1+2/(1+Sqrt))]
EllipticF[I ArcSinh[Sqrt[2/(1+Sqrt)]],(1+Sqrt)/(1-Sqrt)]+I
Sqrt[2/(-1+Sqrt)] EllipticK[(2 Sqrt)/(-1+Sqrt)]

13. The numerical value of this constant:

In:= N[C/.c2rule]

Out= 0.941458+6.66134*10^-16 I

14.  Finally, we define our equation which we want to reverse:

In:=
Clear[finaleq];
finaleq[x_]:=(intfull/.{a->1,C->1})-(C/.c2rule)==-x

In:= finaleq[x]

Out= Sqrt[1/2 (1+Sqrt) (-1-2/(1-Sqrt)) (1+2/(1+Sqrt))]
EllipticF[I ArcSinh[Sqrt[2/(1+Sqrt)]],(1+Sqrt)/(1-Sqrt)]+(I
Sqrt[1/2 (1+Sqrt)] Sqrt[1+(2 y^2)/(1-Sqrt)] Sqrt[1+(2
y^2)/(1+Sqrt)] EllipticF[I ArcSinh[Sqrt[2/(1+Sqrt)]
y],(1+Sqrt)/(1-Sqrt)])/Sqrt[-1+y^2+y^4]==-x

15. Here is the function which solves this equation numerically, to get y =
y[x] (as we only have x = x[y]):

Clear[sol];
sol[x_?NumericQ] := Re[y /. FindRoot[Evaluate[finaleq[x]], {y, 1.5}]]

16. To compare it with something, let us solve the equation numerically by
NDSolve:

In:=
Clear[y];
res = NDSolve[{y''[x]==(1+y'[x]^2+ a
y[x]^4)/y[x]/.a->1,y==1,y'==1},y,{x,0,1}]

During evaluation of In:= NDSolve::ndsz: At x == 0.941458301631319`,
step size is effectively zero; singularity or stiff system suspected. >>

Out= {{y->InterpolatingFunction[{{0.,0.941458}},<>]}}

We see again the same constant value as our C in the analytical solution
as a value where something interesting is happening (a singularity). Here we
define a function representing a solution:

In:= fn = First[y/.res]

Out= InterpolatingFunction[{{0.,0.941458}},<>]

17.  Now we can plot together our two solutions:

Plot[{sol[x], fn[x]}, {x, 0, 1}]

(* Output suppressed *)

Plot[sol[x] - fn[x], {x, 0, 1}]

(* Output suppressed *)

You can see that these are the same curves.

This approach can be generalized to different values of <a> and initial
conditions. While
I did not do the detailed analysis for your conditions, it is quite clear
that there is no simple solution
there because we have an equation like

C-x == h[y,C]

(by h I denoted the integral over y), and since y==y, you have

C  == h[y,C]
C-1 == h[y,C1]],

which clearly has no solution since y=y.

Now, it could be that at x = 0 we take one solution (with plus sign), and at
x =1 - another one (or vice versa),
so that you will have

C  == h[y,C]
C-1 == -  h[y,C1]],

This case I did not look into.

Anyways, you can use the above considerations to perform further analysis of

Hope this helps.

Regards,
Leonid

On Fri, Jul 23, 2010 at 3:13 PM, Sam Takoy <sam.takoy at yahoo.com> wrote:

> I would GREATLY appreciate help symbolically solving the following ODE
> system.
>
> DSolve[y''[x] == (1 + y'[x]^2 + a y[x]^4)/y[x], y, x],
>
> with y = y = 1
>
> and "a" is a constant that satisfies
>
> Integrate[R[x]^2, {x, 0, 1}] = 1/(10*a);
>
> I can't quite pull it off with my current Mathematica skill level.
>
> Many many thanks in advance!
>
> PS: The solution is a smooth function that looks like a parabola.
>
> Thanks again!
>
> Sam
>
>

• Prev by Date: Re: Linking my fortran subroutines with Mathematica
• Next by Date: SetAttributes for entire package
• Previous by thread: Re: Help with an ODE
• Next by thread: Inverse of a six-parameter asymmetric sigmoidal curve