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

```

• Prev by Date: Re: algebraic numbers
• Next by Date: Re: Financial Data - Currencies
• Previous by thread: Re: Difficulty with NDSolve (and DSolve)
• Next by thread: from waveform to third-octave or fft