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}