What's wrong with SetDelayed? The body of the function works alone, while setdelayed to it as $f[x\_]:=func$ gives the wrong answer

So I was writing a function to get the poincare map of a given dynamic system, for example, the duffing oscillator.

I obtained the time series in advance, then take the poincare map based on this time series.

(*define the vector field*)
Duffing[f_, w_] := {v, -0.5*v - x^3 + x + f*Cos[z], w};

(*Obtain the orbits (numerical iteration)*)
OrbitFlow[flow_, vars_, x0_, {t0_, t1_}, opts___] :=
Module[{rules, eq, l = Length[flow], tf, rhs},
rules = Table[vars[[i]] -> vars[[i]][t], {i, 1, l}];
rhs = flow /. rules;
eq = Join[Table[D[vars[[i]] /. rules, t] == rhs[[i]], {i, 1, l}],
Table[vars[[i]] == x0[[i]], {i, 1, l}]];
tf = Table[vars[[i]][t], {i, 1, l}];
sol = NDSolve[eq, tf, {t, t0, t1}, opts]]
(*solving the equation*)

T=2 \[Pi]
sol = OrbitFlow[Duffing[0.39, 1], {x, v, z}, {0, 0, 0}, {0, 1000*T}];

(*to get the set by poincare map*)
poincare[timeseries_,w_]:=Module[{T=2\[Pi]/w,xsol,vsol},
xsol[t_]=x[t]/.timeseries;
vsol[t_]=v[t]/.timeseries;
((Join[xsol[#],vsol[#]]//Evaluate)&) /@Range[800*T,1000*T,T]
];

it doesn't evaluate, and returns:

In:= poincare[sol, 1]
Out:= {{x[1600 \[Pi]], v[1600 \[Pi]]},...}

instead, if I run the code without function-defined 'poincare', it gives evaluated answers,

In:= Module[{T = 2 \[Pi]/1, xsol, vsol}, xsol[t_] = x[t] /. sol;
vsol[t_] = v[t] /. sol;
((Join[xsol[#], vsol[#]] // Evaluate) &) /@ Range[800*T, 1000*T, T]]
Out:= {{-0.508376, 0.294264},...}

Though I know there are plenty of ways to do the same thing, I just wondered why this is wrong?

so it might be something wrong with my settings in NDSolve, in my previous code, I used

NDSolve[eq,{x[t],v[t],z[t]},{~}]

so the output is a list of rules, like:

{{x[t] -> InterpolatingFunction[~],
v[t] -> InterpolatingFunction[~],
z[t] -> InterpolatingFunction[~],
}}

so that evaluation at a specific point may not work, although set another function may solve this but yet cause other problems,

In:= x[0.2]/.sol
Out:= {x[0.2]}

In:= temp[t_] = x[t] /. sol; temp[0.2]
Out:= {0.00754592}

Change the variables in NDSolve may solve this, as

NDSolve[eq,{x,v,z},{~}]

then the solution can be evaluated at some point directly

In:= x[0.2]/.sol
Out:= {0.00754592}

But still, I don't know why this could happen that when setdelayed to it, it will not work.

• Your definition of poincare has one argument and your example uses two arguments. Something is wrong or missing in the code you posted. Nov 29 '21 at 2:55
• yes and thank you, the "1" in poincare[sol,1] was a misstype! Nov 29 '21 at 3:07
• @SuGrape This post is not complete. Could you show one working example of you code to understand your problem? Nov 29 '21 at 3:47
• @AlexTrounev, I posted all the codes, it is based on this article Nov 29 '21 at 4:42

This is actually nontrivial: it has to do with how Mathematica manages scope after substituting in values that match a pattern (as in applying :=).

By inspecting poincare[sol, 1] // TracePrint, we see that after substitution, Mathematica actually uses an expression of the following form:

Module[{T$=2 \[Pi],xsol$,vsol$}, xsol$[t$_] = x[t$] /. {{x[t]->InterpolatingFunction[<<*>>][t], v[t]->InterpolatingFunction[<<*>>][t], z[t]->InterpolatingFunction[<<*>>][t]}};
vsol$[t$_] = v[t$] /. {{x[t]->InterpolatingFunction[<<*>>][t], v[t]->InterpolatingFunction[<<*>>][t], z[t]->InterpolatingFunction[<<*>>][t]}}; (Evaluate[Join[xsol$[#1],vsol$[#1]]] &) /@ Range[800 T$,1000 T$,T$]]

Mathematica has tried to append \$ to the symbol names in an attempt to avoid the very coincidence that you were counting on!

The easiest solution for your case is to extract the InterpolatingFunctions directly as the Head of InterpolatingFunction[<<*>>][t], and set those to xsol and vsol after composing with List (to recover the singleton list structure).

(I'm also not sure why you have Evaluate there, so I took it out, and it seems to work.)

poincare[timeseries_, w_] :=
Module[{T = 2 \[Pi]/w, xsol, vsol},