NDSolve inside Module with local variables

Question

I want to use NDSolve within a Module and output the solution at some times. (Mathematica 10)

All works fine with global variables (symbols)

p[times_List, {k1_, k2_, k3_, k4_}] := Module[{eq},
 eq = {
  m'[t] == -(k2 + k3) m[t] + k1 i[t] + k4 p[t],
  i'[t] == k2 m[t] - k1 i[t],
  p'[t] == k3 m[t] - k4 p[t],
  m[0] == 1, i[0] == 0, p[0] == 0};
 sol = NDSolve[eq, {i[t], m[t], p[t]}, {t, 0, tlist[[-1]]}][[1]];
 pop[t_] = {i[t], m[t], p[t]} /. sol;
 pop /@ times
]

Here an example:

p[Range[0, 5], {1, 1, 2, 3}]

The output is:

{{0., 1., 0.}, {0.312093, 0.404352, 0.283556}, {0.360173, 0.381508, 
 0.25832}, {0.371482, 0.376543, 0.251976}, {0.374165, 0.375366, 
 0.250469}, {0.374802, 0.375087, 0.250111}}

But with local variables (symbols)

p[[times_List, {k1_, k2_, k3_, k4_}] := Module[{eq, t, pop, m, i, p, sol},
 eq = {
  m'[t] == -(k2 + k3) m[t] + k1 i[t] + k4 p[t],
  i'[t] == k2 m[t] - k1 i[t],
  p'[t] == k3 m[t] - k4 p[t],
  m[0] == 1, i[0] == 0, p[0] == 0};
 sol = NDSolve[eq, {i[t], m[t], p[t]}, {t, 0, tlist[[-1]]}][[1]];
 pop[t_] = {i[t], m[t], p[t]} /. sol;
 pop /@ times
]

it does not work anymore:

p[Range[0, 5], {1, 1, 2, 3}]

The output is:

{{i$18625[0], m$18625[0], p$18625[0]}, {i$18625[1], m$18625[1], 
    p$18625[1]}, {i$18625[2], m$18625[2], p$18625[2]}, {i$18625[3], 
    m$18625[3], p$18625[3]}, {i$18625[4], m$18625[4], 
    p$18625[4]}, {i$18625[5], m$18625[5], p$18625[5]}}

Any ideas why this is not working like that? Thank you.

Answer 1

Source of the problem. The problem arises because the symbol t is localized two different times in two different ways, once in Module and once in the definition of pop. Below is a simplified example, with only a single function being solved for, but the same thing happens for OP's system. I print out the form of the solution sol and the definition of pop to show the mismatch of the two symbols represented by t in the code. The mismatch means that the solution is not substituted for x[t] in the definition of pop. Because the body of the definition of pop depends on Module variables, the function parameter t is remaned t$. As it says in the tutorial, Variables in Pure Functions and Rules, Mathematica does this more liberally than is strictly necessary. Because of this, we have x$311828[t$] in the body of pop and x$311828[t$311828] in rule in sol, and they do not match.

p[times_List, {k1_, k2_, k3_, k4_}] := 
 Module[{eq, t, pop, x, sol}, 
  sol = NDSolve[{x'[t] == x[t], x[0] == 1}, {x[t]}, {t, 0, 5}][[1]];
  Print[sol];
  pop[t_] = {x[t]} /. sol;
  Print[Definition[pop]];
  pop /@ times]

p[Range[0, 5], {1, 1, 2, 3}]

(*
  {x$311828[t$311828] -> InterpolatingFunction[{{0., 5.}}, <>][t$311828]}
*)
(*
  Attributes[pop$311828]={Temporary}

  pop$311828[t$_]={x$311828[t$]}
*)
(*
  {{x$311828[0]}, {x$311828[1]}, {x$311828[2]}, {x$311828[3]}, {x$311828[4]}, {x$311828[5]}}
*)

Personally, I have come to prefer to call NDSolve in the form

NDSolve[eqns, x, {t, 0, 5}]

with x instead of x[t], and this is a way to solve the OP's problem. NDSolveValue is an attractive alternative, but I have come to prefer to keep the symbolic expressions and equations separate from the approximate solutions. The call above produces a rule of the form

  {x -> InterpolatingFunction[{{0., 5.}}, <>]}

which replaces a head x in x[t] or x[1.2] by a function, which will automatically evaluate if the argument is a number. This means that it does not matter how t is localized; the substitution will work.

Changing one line fixes the code. Replace the NDSolve line with this:

 sol = NDSolve[eq, {i, m, p}, {t, 0, tlist[[-1]]}][[1]];

Further remarks. There are several advantages to using NDSolve this way. The main disadvantage is writing /. sol a lot.

One advantage is that the symbol x (or i, m, and p) remains undefined and can be used for symbolic manipulation and analysis.

The substitution /. sol works on derivatives, x'[t] etc., which is awkward to accomplish if the rule is in the form x[t] -> ....

You can check the residuals to see how good the solution is. For instance,

eq /. Equal -> Subtract /. sol /. t -> RandomReal[{0, 5}, 10]

or one can use Table or NIntegrate instead of random sampling.

Plotting is just as easy, and if the integration stops early and you don't know the domain exactly, you can get it out of sol:

Plot[Evaluate[x[t] /. sol], Evaluate[Flatten[{t, x["Domain"] /. sol}]]