Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Difficulty with NDSolve (and DSolve)

• To: mathgroup at smc.vnet.net
• Subject: [mg106167] Re: Difficulty with NDSolve (and DSolve)
• From: Mark McClure <mcmcclur at unca.edu>
• Date: Sun, 3 Jan 2010 03:43:31 -0500 (EST)
• References: <201001021004.FAA07409@smc.vnet.net>

On Sat, Jan 2, 2010 at 5:04 AM, KK <kknatarajan at 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

• References:
  • Difficulty with NDSolve (and DSolve)
    • From: KK <kknatarajan@yahoo.com>