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>
- Problem with a 1st order IV ODE (nonlinear)