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

Re: NDSolve does not work in shooting method

• To: mathgroup at smc.vnet.net
• Subject: [mg44159] Re: NDSolve does not work in shooting method
• From: Rob Knapp <rknapp at wolfram.com>
• Date: Fri, 24 Oct 2003 04:24:46 -0400 (EDT)
• References: <bn8du7$nfi$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Alexander Koulidis wrote:
> This comes from a newbie in Mathematica, so have some mercy....
> 
> I'm trying to reproduce some examples of the shooting method that have been
> published in this forum, but Mathematica 5.0 shows some strange behavior
> with NDSolve.
> 
> In http://forums.wolfram.com/mathgroup/archive/2001/Aug/msg00128.html Brian
> wrote:
> 
> The following code implements the shooting technique
> 
> system[a_] := {y1'[x] == y2[x], y2'[x] == -Sin[y2[x]] + Cos[5x],
> y1[-1] == 0,
>       y2[-1] == a};
> myODEsoln[a_] := NDSolve[system[a], {y1[x], y2[x]}, {x, -1, 1}]
> yend[a_] := (y1[x] /. myODEsoln[a]) /. x -> 1
> bc = FindRoot[First[yend[a]] == 0, {a, -2, 2}];
> Plot[Evaluate[y1[x] /. myODEsoln[a /. bc]], {x, -1, 1},
>     AxesLabel -> {"x", "y1(x)"}];
> 
> 
> While the code was working fine in v. 4.2, now it produces warnings of
> "NDSolve::ndinnt :" or "RuleDelayed::rhs".
> Can someone test this code in v5.0? Does anyone have some idea of what is
> happening and should I implement a typical shooting method?
> Thanks in advance.
> Alexander Koulidis, Athens, Greece
> 
> 

The strange behavior is not in NDSolve, but arises from the way 
expressions are evaluated, in particular in FindRoot.  As of Mathematica 
5, what FindRoot (also FindMinimum, NIntegrate, NSum) does is give the 
variables local definitions, then evaluates the function (or equations 
for FindRoot) with the variables having no values: i.e. symbolically.


This means that in your example, the equation FindRoot is working with 
is not that you were expecting:

In[1]:= system[a_] := {y1'[x] == y2[x], y2'[x] == -Sin[y2[x]] + 
Cos[5x],y1[-1] == 0,      y2[-1] == a};

In[2]:= myODEsoln[a_] := NDSolve[system[a], {y1[x], y2[x]}, {x, -1, 1}];

In[3]:= yend[a_] := (y1[x] /. myODEsoln[a]) /. x -> 1;

In[4]:= yend[a]

.... various messages elided ...


Out[4]= y1[1] /.

 >    NDSolve[{y1'[1] == y2[1], y2'[1] == Cos[5] - Sin[y2[1]], y1[-1] == 0,

 >      y2[-1] == a}, {y1[1], y2[1]}, {1, -1, 1}]

In[5]:= First[%]

Out[5]= y1[1]

.... so the equation FindRoot tries to work with is

y1[1] == 0, which is not really what you want to solve.


The real problem here is that the function yend does not do what you 
intend unless the argument a is a number.  If you simply attach this 
restriction to the definition to delay evaluation until it is 
appropriate to evaluate, then it works as you expect:

In[6]:= Clear[yend]

In[7]:= yend[a_?NumberQ] := First[(y1[x] /. myODEsoln[a]) /. x -> 1]

Now yend does not evaluate with symbolic a

In[8]:= yend[a]

Out[8]= yend[a]

(Note I included the First in the definition of yend, since 
First[yend[a]] would evaluate to a with non number values of a)

It does evaluate with numbers.

In[9]:= yend[1]

Out[9]= 0.715203

.... and FindRoot works just fine with this definition.

In[10]:= FindRoot[yend[a] == 0, {a, -2,2}]

Out[10]= {a -> 0.212199}