2
$\begingroup$

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.

(All the code is based on this article)

(*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]][0] == 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.


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

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}, 
  xsol = List@*Head[First[x[t] /. timeseries]]; 
  vsol = List@*Head[First[v[t] /. timeseries]];
  (Join[xsol[#], vsol[#]]) & /@ Range[800*T, 1000*T, T]]

(There are other more in-depth ways to avoid those scoping conflicts, but given that this is not essential to what you're doing, I think it makes more sense to just go directly for what you want! But let me know if you're interested!)

$\endgroup$
1
  • $\begingroup$ thank you sooo much!!! It took me a while to understand, and yes I am interested, it would be very generous of you to show me more references. $\endgroup$
    – Su Grape
    Nov 29 '21 at 7:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.