Prev: calling mathematica functions from c++
Next: A question on plot3D
From: Sam Takoy on 23 Jul 2010 07:13
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[0] = y[1] = 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

From: Sam Takoy on 26 Jul 2010 07:11
On 7/23/2010 7:13 AM, Sam Takoy 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[0] = y[1] = 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
>

Hi,

May I just add that I would be very happy with a solution to

y''[x] == (1 + y'[x]^2 + y[x]^4)/y[x]

In other words, forget "a". I believe that the case a=1 has full generality.

Once again, I would greatly appreciate and help with solving this system.

Many many thanks in advance!

Sam

From: Leonid Shifrin on 27 Jul 2010 04:15
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[0]=1, y'[0]=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[77]:= eq = y''[x]==(1+y'[x]^2+ a y[x]^4)/y[x]

Out[77]= (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[81]:= eq1=eq/.{y'[x]->f[y],y''[x]:> 1/2D[f[y]^2,y],y[x]-> y}

Out[81]= 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[88]:= eq2 =Apart/@(eq1/.{f[y]
(f^\[Prime])[y]->1/2D[g[y],y],f[y]^2->g[y]})

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

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

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

Out[91]= -1 + a y^4 + y^2 C[1]

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

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

Out[93]= True

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

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

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

(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[97]:= HoldForm[-x+C[2]==-Integrate[Sqrt[a t^4+C[1]t^2-1],{t,y,Infinity}]]
Out[97]= (* 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[107]:= int = Integrate[1/Sqrt[a t^4+C[1] t^2-1],t]

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

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

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

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


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

In[110]:= inty = int /. t -> y

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

9. The final result for the integral then:

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


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

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

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

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


11. We have then :

In[119]:= -x+C[2]==intfull

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

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

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

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

13. The numerical value of this constant:

In[118]:= N[C[2]/.c2rule]

Out[118]= 0.941458+6.66134*10^-16 I

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

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

In[133]:= finaleq[x]

Out[133]= Sqrt[1/2 (1+Sqrt[5]) (-1-2/(1-Sqrt[5])) (1+2/(1+Sqrt[5]))]
EllipticF[I ArcSinh[Sqrt[2/(1+Sqrt[5])]],(1+Sqrt[5])/(1-Sqrt[5])]+(I
Sqrt[1/2 (1+Sqrt[5])] Sqrt[1+(2 y^2)/(1-Sqrt[5])] Sqrt[1+(2
y^2)/(1+Sqrt[5])] EllipticF[I ArcSinh[Sqrt[2/(1+Sqrt[5])]
y],(1+Sqrt[5])/(1-Sqrt[5])])/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[139]:=
Clear[y];
res = NDSolve[{y''[x]==(1+y'[x]^2+ a
y[x]^4)/y[x]/.a->1,y[0]==1,y'[0]==1},y,{x,0,1}]

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

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

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

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

Out[141]= 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[2]-x == h[y,C[1]]

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

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

which clearly has no solution since y[0]=y[1].

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[2] == h[y[0],C[1]]
C[2]-1 == - h[y[1],C1]],

This case I did not look into.

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

Hope this helps.

Regards,
Leonid




On Fri, Jul 23, 2010 at 3:13 PM, Sam Takoy <sam.takoy(a)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[0] = y[1] = 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
>
>
 | 
Pages: 1
Prev: calling mathematica functions from c++
Next: A question on plot3D