Re: Difficulty with NDSolve (and DSolve)
- To: mathgroup at smc.vnet.net
- Subject: [mg106157] Re: Difficulty with NDSolve (and DSolve)
- From: danl at wolfram.com
- Date: Sun, 3 Jan 2010 03:41:36 -0500 (EST)
- References: <201001021004.FAA07409@smc.vnet.net>
> Hi,
>
> I am trying to numerically solve a differential equation. However, I
> am encountering difficulty with Mathematica in generating valid
> numerical solutions even for a special case of that equation.
>
> The differential equation for the special case is:
> F'[x] == - (2-F[x])^2/(1-2 x + x F[x]) and
> F[1]==1.
> These equations are defined for x in (0,1). Moreover, for my context,
> I am only interested in solutions with F[x] in the range <1.
>
> Even before I used Mathematica, I had computed the solution to the
> differential equation as the solution to the following : x (2-F[x]) -
> Log[2-F[x]]==1. For any given x in (0,1), there are two values of F
> [x] that satifsy the equation. One of them is always less than 1, and
> the other is F[x]= 2 + (ProductLog[E^(I (I + \[Pi])) x])/x which is
> always greater than 1. For example, when x=0.85, the solutions are F
> [x]=0.2979 and F[x]=1.32407. As mentioned earlier, I am only
> interested in the first solution.
>
> Both DSolve and NDSolve seem to provide only the second solution and
> not the first one. When using DSolve, it gives out a warning that
> there may be multiple solutions but that's about it. Is there an
> option or something that I can set with NDSolve (or even DSolve) to
> generate the solutions of interest to me?
>
> Below is the relevant set of Mathematica code and the corresponding
> outputs. I would greatly appreciate your help in this regard.
>
> Thanks,
> KK
>
>
> ------------------
> In[1]:= DSolve[{ -((-2 + F[x])^2/(1 - 2 x + x F[x])) == F'[x],
> F[1] == 1}, F[x], x]
>
> During evaluation of In[1]:= InverseFunction::ifun: Inverse functions
> are being used. Values may be lost for multivalued inverses. >>
>
> During evaluation of In[1]:= Solve::ifun: Inverse functions are being
> used by Solve, so some solutions may not be found; use Reduce for
> complete solution information. >>
>
> During evaluation of In[1]:= Solve::ifun: Inverse functions are being
> used by Solve, so some solutions may not be found; use Reduce for
> complete solution information. >>
>
> Out[1]= {{F[x] -> (2 x + ProductLog[E^(I (I + \[Pi])) x])/x}}
> ------------------
>
> (*Sample solution*)
> In[2]:= Solve[(1 == (x (2 - F[x]) - Log[2 - F[x]] /. F[x] -> F) /.
> x -> 0.9), F]
>
> During evaluation of In[2]:= InverseFunction::ifun: Inverse functions
> are being used. Values may be lost for multivalued inverses. >>
>
> During evaluation of In[2]:= Solve::ifun: Inverse functions are being
> used by Solve, so some solutions may not be found; use Reduce for
> complete solution information. >>
>
> Out[2]= {{F -> 0.297987}, {F -> 1.32407}}
>
> --------------------
> (*Please note I am using F[0.9999]=1 as the initial condition in the
> NDSolve below. \
> Otherwise, I end up with a warning that "Infinite expression 1/0. \
> encountered. I had also used the same initial condition with DSolve \
> and the results were similar.*)
>
> In[3]:= NumSol =
> NDSolve[{ -((-2 + F[x])^2/(1 - 2 x + x F[x])) == F'[x], F[.9999] ==
> 1}, F, {x, 0, 1}]
>
> During evaluation of In[3]:= NDSolve::ndsz: At x ==
> 0.9999000049911249`, step size is effectively zero; singularity or
> stiff system suspected. >>
>
> Out[156]= {{F -> \!\(\*
> TagBox[
> RowBox[{"InterpolatingFunction", "[",
> RowBox[{
> RowBox[{"{",
> RowBox[{"{",
> RowBox[{"0.`", ",", "0.9999000049911249`"}], "}"}], "}"}],
> ",", "\<\"<>\"\>"}], "]"}],
> False,
> Editable->False]\)}}
> ---
> (* I have omitted the plot here but it generates only the second
> result*)
> In[157]:= Plot[{Evaluate[F[x] /. NumSol]}, {x, .33, 1}]
> -----------------------
>
You can forge the DSolve result by explicitly changing branches of
ProductLog.
sol = DSolve[{-((-2 + F[x])^2/(1 - 2 x + x F[x])) == F'[x],
F[1] == 1}, F[x], x];
ff = F[x] /. sol[[1]] /. ProductLog[aa_] :> ProductLog[-1, aa];
Evaluating at .8 (not .85) gives the stated result you have in mind.
In[75]:= Chop[ff2 /. x -> .9]
Out[75]= 0.29798710178932075
For NDSolve, the issue is to tilt the initial value in the other direction.
In[80]:= -((-2 + y)^2/(1 - 2 x + x y)) /. {y -> 1, x -> .9999}
Out[80]= -10000.0000000011
You want a result that is increasing as you approach from the left, so it
requires a positive derivative there.
In[81]:= -((-2 + y)^2/(1 - 2 x + x y)) /. {x -> 1, y -> .9999}
Out[81]= 10002.0001000011
So instead of setting F[.9999] to 1, you might set F[1] to .9999.
In[82]:= NumSol =
NDSolve[{-((-2 + F[x])^2/(1 - 2 x + x F[x])) == F'[x],
F[1] == .9999}, F, {x, 0, 1}];
During evaluation of In[82]:= NDSolve::mxst: Maximum number of 10000 steps
reached at the point x == 1.2113589086340555`*^-207. >>
In[83]:= ff = F /. NumSol[[1]];
In[84]:= ff[.9]
Out[84]= 0.29798716075095194
Daniel Lichtblau
Wolfram Research
- References:
- Difficulty with NDSolve (and DSolve)
- From: KK <kknatarajan@yahoo.com>
- Difficulty with NDSolve (and DSolve)