Re: Problem with a 1st order IV ODE (nonlinear)

*To*: mathgroup at smc.vnet.net*Subject*: [mg102581] Re: [mg102564] Problem with a 1st order IV ODE (nonlinear)*From*: Leonid Shifrin <lshifr at gmail.com>*Date*: Sat, 15 Aug 2009 05:34:03 -0400 (EDT)*References*: <200908140959.FAA01407@smc.vnet.net>

Hi Virgil, while I don't have an authorative explanation of the behavior of NDSolve in your case, I do have an alternative approach to suggest. First of all, I have to disappoint you: the solutions you (or, more precisley, DSolve) obtained look nice, but don't seem to work for me: In[1] = Clear[R, k]; sol = DSolve[{h'[t] == 1/(h[t] (2 R - h[t])) - k, h[0] == 0} /. {R -> 10, k -> 0.01}, h[t], t] // FullSimplify Out[1] = {{h[t]->-0.005 t-0.005 Sqrt[t (4000.+t)]},{h[t]->-0.005 t+0.005 Sqrt[t (4000.+t)]}} In[2] = fn1 = FullSimplify[D[#, t] - 1/(# (2 R - #) - k) /. {R -> 10., k -> 0.01}] &[(h[ t] /. sol)[[1]]] Out[2] = -0.005+(-10.-0.005 t)/Sqrt[t (4000.+t)]+1/(0.01+0.1 Sqrt[t (4000.+t)]+t (0.2+0.00005 t+0.00005 Sqrt[t (4000.+t)])) In[3] = fn2 = FullSimplify[D[#, t] - 1/(# (2 R - #) - k) /. {R -> 10., k -> 0.01}] &[(h[t] /. sol)[[2]]] Out[3] = -0.005+(10.+0.005 t)/Sqrt[t (4000.+t)]+1/(0.01-0.1 Sqrt[t (4000.+t)]+t (0.2+0.00005 t-0.00005 Sqrt[t (4000.+t)])) You can plot these functions: In[4] = Plot[{fn1, fn2}, {t, 0, 1}] and see that they are non-zero, although small. Now, I don't know why DSolve does this, but we can simply integrate your equaion: In[5] = Clear[R, k]; Integrate[(h (2 R - h)) - k, h] Out[5] = -(h^3/3) - h k + h^2 R The resulting algebraic equation t == -(h^3/3) - h k + h^2 R must be solved for h. We solve now: In[6] = sols = Solve[t + c == -(h^3/3) - h k + h^2 R, h] Out[6] = {{h->R-(2^(1/3) (9 k-9 R^2))/(3 (-81 c-81 k R+54 R^3+Sqrt[4 (9 k-9 R^2)^3+(-81 c-81 k R+54 R^3-81 t)^2]-81 t)^(1/3))+(-81 c-81 k R+54 R^3+Sqrt[4 (9 k-9 R^2)^3+(-81 c-81 k R+54 R^3-81 t)^2]-81 t)^(1/3)/(3 2^(1/3))}, {h->R+((1+I Sqrt[3]) (9 k-9 R^2))/(3 2^(2/3) (-81 c-81 k R+54 R^3+Sqrt[4 (9 k-9 R^2)^3+(-81 c-81 k R+54 R^3-81 t)^2]-81 t)^(1/3))-((1-I Sqrt[3]) (-81 c-81 k R+54 R^3+Sqrt[4 (9 k-9 R^2)^3+(-81 c-81 k R+54 R^3-81 t)^2]-81 t)^(1/3))/(6 2^(1/3))}, {h->R+((1-I Sqrt[3]) (9 k-9 R^2))/(3 2^(2/3) (-81 c-81 k R+54 R^3+Sqrt[4 (9 k-9 R^2)^3+(-81 c-81 k R+54 R^3-81 t)^2]-81 t)^(1/3))-((1+I Sqrt[3]) (-81 c-81 k R+54 R^3+Sqrt[4 (9 k-9 R^2)^3+(-81 c-81 k R+54 R^3-81 t)^2]-81 t)^(1/3))/(6 2^(1/3))}} We can check back that all solutions satisfy the original diff. equation: In[7] = FullSimplify[D[#, t] - 1/((# (2 R - #)) - k) &[h /. #]] & /@ sols Out[7] = {0, 0, 0} Neither of the solutions, however, satisfies the initial condition for all values of the parameters R and t - rather, for different {R,t} values you pick those that satisfy. Namely, we have to solve for the integration constant <c>, and I'm afraid, numerically, for your particular values of parameters R,k: In[8] = construles = NSolve[((h /. #) /. {t -> 0, R -> 10., k -> 0.01}) == 0, c] & /@ sols Out[8] = {{},{{c->3.40612*10^-14-1.2179*10^-17 I}},{{c->3.40612*10^-14+1.2179*10^-17 I}}} The result suggests that for the first solution, there is no solution for <c> (for these particular values of the parameters), while for the other two it is effectively zero (or so it seems, see below!) Here are your two solutions for your particular parameters then: In[9] = {ff1, ff2} = h /. Rest[sols] /. {c -> 0, R -> 10, k -> 0.01}\ Out[9] = { 10-(188.969+327.304 I)/(53991.9+Sqrt[-2.91513*10^9+(53991.9-81 t)^2]-81 t)^(1/3)-((1-I Sqrt[3]) (53991.9+Sqrt[-2.91513*10^9+(53991.9-81 t)^2]-81 t)^(1/3))/(6 2^(1/3)), 10-(188.969-327.304 I)/(53991.9+Sqrt[-2.91513*10^9+(53991.9-81 t)^2]-81 t)^(1/3)-((1+I Sqrt[3]) (53991.9+Sqrt[-2.91513*10^9+(53991.9-81 t)^2]-81 t)^(1/3))/(6 2^(1/3))} The first one satisfies the initial condition well: In[10] = ff1 /. t -> 0 Out[10] = 5.14877*10^-12-5.32907*10^-15 I The second one does not In[11] = ff2 /. t -> 0 Out[11] = 0.00100003+5.32907*10^-15 I The reason is that the equation is singular, and we in fact should not have ignored the small integration constant in this case. Defining ff2Correct as In[12] = ff2Correct = h /. sols[[3]] /. {c -> (c /. construles[[3, 1]]), R -> 10, k -> 0.01} Out[12] = 10-(188.969-327.304 I)/((53991.9-9.86503*10^-16 I)+Sqrt[-2.91513*10^9+((53991.9-9.86503*10^-16 I)-81 t)^2]-81 t)^(1/3)-((1+I Sqrt[3]) ((53991.9-9.86503*10^-16 I)+Sqrt[-2.91513*10^9+((53991.9-9.86503*10^-16 I)-81 t)^2]-81 t)^(1/3))/(6 2^(1/3)) we have In[13] = ff2Correct /. t -> 0 Out[13] = -5.75451*10^-12+3.55271*10^-15 I You can check that solutions satisfy the original differential equation In[14] = checks = D[#, t] - 1/(# (2 R - #) - k) /. {R -> 10., k -> 0.01} & /@ {ff1, ff2Correct}; In[15] = Plot[Re[checks], {t, 0, 1}] (Output omitted) In[16] = Plot[Im[checks], {t, 0, 1}] (Output omitted) It is interesting that the solutions we obtained are pretty close to the one you got from DSolve, and seem to differ from them in a term linear in t. You can use the above procedure for other values of R,k. In general, also for the first solution the value of the integration constant should not be neglected, even if very small. Also, for different values of R,k, I expect different combinations of the above solutions to satisfy initial condition. Hope this helps. Regards, Leonid On Fri, Aug 14, 2009 at 1:59 PM, Virgil Stokes <vs at it.uu.se> wrote: > I am using Mathematica 7.0 on a Win2K platform and noticed that when I > execute the following: > > R = 10; > k = 0.01; > sol = DSolve[{h'[t] == 1/(h[t] (2 R - h[t])) - k, h[0] == 0}, h[t], t] > // FullSimplify > > I get two possible solutions: > {{h[t] -> -0.005 t - 0.005 Sqrt[t (4000. + t)]}, {h[t] -> -0.005 t + > 0.005 Sqrt[t (4000. + t)]}} > > which, I believe are correct. However, if I try to get an analytical > solution in terms of R and k, > > Clear[R, k] > sol = DSolve[{h'[t] == 1/(h[t] (2 R - h[t])) - k, h[0] == 0}, h[t], t] > // FullSimplify > > I get the following two output messages: > Solve::tdep: The equations appear to involve the variables to be solved > for in an essentially non-algebraic way. > DSolve::bvnul: For some branches of the general solution, the given > boundary conditions lead to an empty solution. > > Note, that R > 0, and k >= 0.. > Is there anyway that I can get an analytical solution to this problem > for these conditions? > > --V. Stokes > > >

**References**:**Problem with a 1st order IV ODE (nonlinear)***From:*Virgil Stokes <vs@it.uu.se>