Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Wrapping NDSolve within a function

• To: mathgroup at smc.vnet.net
• Subject: [mg126355] Re: Wrapping NDSolve within a function
• From: Patrick Scheibe <pscheibe at trm.uni-leipzig.de>
• Date: Thu, 3 May 2012 22:21:45 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <201205020944.FAA19622@smc.vnet.net>

Hi,

this is because of the lexical scoping which is done by Module. The g, p
and q are renamed to not clash with some globaly defined variables. You
could prevent this by using Block which uses dynamic scoping and does
not rename the variables.

predState[rg_, rp_, rh_, Mg_, Mp_, Mh_] := 
 Module[{parmList3woC, parmListND3woC, A, B, Rb, Rc, solND}, 
  parmListND3woC = {A -> Mg/rg, B -> Mp/rg, C -> Mh/rg, Rb -> rp/rg, 
    Rc -> rh/rg};
  Block[{g, p, q}, 
   NDSolve[{g'[t] == g[t]*(1 - g[t] - p[t] - q[t]) - A*g[t]*g[t] /. 
      parmListND3woC, 
     p'[t] == Rb*p[t]*(1 - p[t] - q[t]) - B*p[t]*g[t] /. 
      parmListND3woC, 
     q'[t] == Rc*q[t]*(1 - q[t]) - C*q[t]*g[t] /. parmListND3woC, 
     g[0] == 0.5, p[0] == 0.5, q[0] == 0.5}, {g, p, q}, {t, 0, 1000}]
   ]
  ]

But be aware what happens if you set for instance g=3; before calling
your function!

Why don't you return a function? With this you could just *use* your
InterpolatingFunctions:

predState[rg_, rp_, rh_, Mg_, Mp_, Mh_] := 
 Module[{parmList3woC, parmListND3woC, A, B, Rb, Rc, g, p, q}, 
  parmListND3woC = {A -> Mg/rg, B -> Mp/rg, C -> Mh/rg, Rb -> rp/rg, 
    Rc -> rh/rg};
  With[{sol = {g, p, q} /. 
      First@NDSolve[{g'[t] == 
           g[t]*(1 - g[t] - p[t] - q[t]) - A*g[t]*g[t] /. 
          parmListND3woC, 
         p'[t] == Rb*p[t]*(1 - p[t] - q[t]) - B*p[t]*g[t] /. 
          parmListND3woC, 
         q'[t] == Rc*q[t]*(1 - q[t]) - C*q[t]*g[t] /. parmListND3woC, 
         g[0] == 0.5, p[0] == 0.5, q[0] == 0.5}, {g, p, q}, {t, 0, 
         1000}]
    },
   Function[{t}, Through[sol[t]]]
   ]
  ]

out = predState[1.5, 0.4, 0.2, 0.1, 0.0002, 0.0099];
Plot[out[t], {t, 0, 10}]

Cheers
Patrick


On Wed, 2012-05-02 at 05:44 -0400, bbeckage wrote:
> When I try to return the interpolating function produced by NDSolve from 
> within a function, the object returned has an unexpected $7360 appended, 
> e.g.,
> 
> out = predState[1.5, 0.4, 0.2, 0.1, 0.0002, 0.0099]
> 
> (where predState is defined further below) results in
> 
> =
> {{g$7360->InterpolatingFunction[{{0.,1000.}},<>],p$7360->InterpolatingFunction[{{0.,1000.}},<>],q$7360->InterpolatingFunction[{{0.,1000.}},<>]}}
> 
> Note the g$7360 rather than just g.  If NDSolve is not wrapped within 
> the function, it returns a plain 'g', i.e., g->InterpolatingFunction.... 
>  The appended $7360 makes it difficult to use the interpolating function 
> as I can't reference it within other functions as the integer changes 
> with each function call, e.g., g$7360[10] /. out, then g$7370[10] /. 
> out, rather than being able to access it using the expected g[10]/.out.
> 
> Why  is this $7360 appended to g?  How can NDSolve be wrapped in a 
> function, but be made to return a plain g?
> 
> Thanks for your help.
> 
> Best wishes,
> Brian
> 
> 
> 
> predState[rg_, rp_, rh_, Mg_, Mp_, Mh_] :=
>  Module[{parmList3woC, parmListND3woC, A, B, Rb, Rc, g, p, q, solND},
>   parmListND3woC = {A -> Mg/rg, B -> Mp/rg, C -> Mh/rg, Rb -> rp/rg,
>     Rc -> rh/rg};
>   solND =
>    NDSolve[{
>      g'[t] == g[t]*(1 - g[t] - p[t] - q[t]) - A*g[t]*g[t] /. parmListND3woC,
>      p'[t] == Rb*p[t]*(1 - p[t] - q[t]) - B*p[t]*g[t] /. parmListND3woC,
>      q'[t] == Rc*q[t]*(1 - q[t]) - C*q[t]*g[t] /. parmListND3woC, g[0] == 0.5,
>       p[0] == 0.5, q[0] == 0.5}, {g, p, q}, {t, 0, 1000}];
>   Return[solND]
>   ]
> 
> 
> 
> 
> 
> 
> 
> 

• References:
  • Wrapping NDSolve within a function
    • From: bbeckage <Brian.Beckage@uvm.edu>