Re: InterpolatingFunction
- To: mathgroup at smc.vnet.net
- Subject: [mg5020] Re: [mg4972] InterpolatingFunction
- From: jpk at apex.mpe.FTA-Berlin.de (Jens-Peer Kuska)
- Date: Sat, 19 Oct 1996 16:40:27 -0400
- Sender: owner-wri-mathgroup at wolfram.com
----------
X-Sun-Data-Type: text
X-Sun-Data-Description: text
X-Sun-Data-Name: text
X-Sun-Content-Lines: 59
X-Sun-Charset: us-ascii
Hi,
I append a package
Poincare.m that does all what You want. How ever it is
a bit slow because is was written for explaning the
numerical technique not for calculating *real* Poincare
sections that have usual some hundred or thousend points.
For that case You can have my sample program henonml thats
compute the Poincare section for the Henon-Heiles system
in a Mathlink program. Contact me directly is You want the
C-Code a cover package to launch the ML-program and a sample
notebook.
Hope that helps
Jens
BeginPackage["Poincare`"]
PoincareSection::usage=
"PoincareSection[traj,plane,{t,t0,tE}] calculates the
Poincare Section of an autonomus dynamical system. The
trajectorie of the system traj is the solution returned
by NDSolve. PoicareSection returns the list of
intersection points with the plane plane. The
time interval t in [t0,tE] is searched for possible
intersection points. A optional fourth argument
can be given to increase th initial number of 500
subdivisions of the time interval."
Begin["`Private`"]
PoincareSection[traj_,plane:Equal[left_,right_],
time:{t_Symbol,t0_?NumberQ,tE_?NumberQ},
subdiv_:500]:=
Module[{sectlst={},dt=N[(tE-t0)/subdiv],
tau=t0,sect,
eqn,pplane=plane /. traj,
interpol= #[[2]] & /@ Flatten[traj],
sign1,sign2},
eqn=(left-right) /. traj;
sign1=Sign[eqn /. t->t0 ];
While[tau<tE,
tau+=dt;
If[tau>tE, tau=tE];
sign2=Sign[eqn /. t->tau];
If[sign1=!=sign2,
sect=FindRoot[pplane,{t,tau-dt,tau}];
If[sect =!=$Failed,
AppendTo[sectlst,interpol /. sect]];
];
sign1=sign2;
];
sectlst
]
End[]
EndPackage[]
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: POINCARE.M
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 42
BeginPackage["Poincare`"]
PoincareSection::usage=
"PoincareSection[traj,plane,{t,t0,tE}] calculates the
Poincare Section of an autonomus dynamical system. The
trajectorie of the system traj is the solution returned
by NDSolve. PoicareSection returns the list of
intersection points with the plane plane. The
time interval t in [t0,tE] is searched for possible
intersection points. A optional fourth argument
can be given to increase th initial number of 500
subdivisions of the time interval."
Begin["`Private`"]
PoincareSection[traj_,plane:Equal[left_,right_],
time:{t_Symbol,t0_?NumberQ,tE_?NumberQ},
subdiv_:500]:=
Module[{sectlst={},dt=N[(tE-t0)/subdiv],
tau=t0,sect,
eqn,pplane=plane /. traj,
interpol= #[[2]] & /@ Flatten[traj],
sign1,sign2},
eqn=(left-right) /. traj;
sign1=Sign[eqn /. t->t0 ];
While[tau<tE,
tau+=dt;
If[tau>tE, tau=tE];
sign2=Sign[eqn /. t->tau];
If[sign1=!=sign2,
sect=FindRoot[pplane,{t,tau-dt,tau}];
If[sect =!=$Failed,
AppendTo[sectlst,interpol /. sect]];
];
sign1=sign2;
];
sectlst
]
End[]
EndPackage[]