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