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>