Difficulty with NDSolve (and DSolve)
in [Mathematica]
|
From: KK on 2 Jan 2010 05:04 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}] -----------------------
From: danl on 3 Jan 2010 03:38 > 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
From: Mark McClure on 3 Jan 2010 03:39 On Sat, Jan 2, 2010 at 5:04 AM, KK <kknatarajan(a)yahoo.com> wrote: > 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
|
Pages: 1 Prev: Financial Data - Currencies Next: NDSolve problem with switching equations |