[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
Prev by Date:
**Re: Financial Data - Currencies**
Next by Date:
**Re: More /.{I->-1} craziness**
Previous by thread:
**Difficulty with NDSolve (and DSolve)**
Next by thread:
**Re: Difficulty with NDSolve (and DSolve)**
| |