[Date Index]
[Thread Index]
[Author Index]
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[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 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[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
>
>
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**
| |