Re: Difficulty with NDSolve (and DSolve)
da...@wolfram.com
>
> 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
Mark McClure
> 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.
Of course, the basic problem is that the solution is not one-to-one in
a neighborhood of your initial condition. Thus, let's set up an
initial condition y0 at x0=1/2 and try to choose y0 so that the
solution passes through (1,1) - kind of like the shooting method for
solving boundary value problems. We can set up a function that tells
us the value of the solution at x=1 as a function of the initial
condition at x0=1/2 like as follows:
In[1]:= valAt1[y0_?NumericQ] := Quiet[Check[(F[x] /.
NDSolve[{F'[x] == -(2 - F[x])^2/(1 - 2 x + x F[x]),
F[1/2] == y0}, F[x], {x, 1/2, 1}][[1, 1]]) /. x -> 1,
bad]];
In[2]:= {valAt1[-4], valAt1[-3]}
Out[2]= {0.208967, bad}
Note that we use Check to check for error messages in case that x=1 is
outside the range of the solution. The computation shows that this
singularity between y0=-4 and y0=-3. We can use a bisection method to
find more precisely where this happens like so:
In[3]:= a = -4;
In[4]:= b = -3;
In[5]:= Do[
c = (a + b)/2;
If[valAt1[c] === bad, b = c, a = c],
{50}];
In[6]:= N[a]
Out[6]= -3.35669
Now, we find the solution taking y0=a.
In[7]:= Clear[F];
In[8]:= F[x_] = F[x] /.
NDSolve[{F'[x] == -(2 - F[x])^2/(1 - 2 x + x F[x]),
F[1/2] == a}, F[x], {x, 0.001, 1}][[1]];
In[9]:= F[1]
Out[9]:= 1.
The computation shows that it worked.
Hope that helps,
Mark McClure