Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.math.mathematica

Topic: Wrapping NDSolve within a function
Replies: 2   Last Post: May 3, 2012 10:22 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Patrick Scheibe

Posts: 327
Registered: 5/22/07
Re: Wrapping NDSolve within a function
Posted: May 3, 2012 10:22 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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]
> ]
>
>
>
>
>
>
>
>







Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.