Master di II livello in Scienza e Tecnologia Spaziale
METODI NUMERICI PER L' ASTRONOMIA

Il software Mathematica per il calcolo simbolico, numerico e grafico.

Tor Vergata, 8 febbraio 2015

Corrado Falcolini

Università degli Studi Roma Tre

(* Questo file si può far girare lanciando Evaluate Notebook dalla barra Evaluation *)

Problema 3: visualizzazione e calcolo di orbite per la Standard Map.

La Standard Map

Un esempio di sistema dinamico discreto nel piano è la cosiddetta Mappa Standard definita sul toro   (T^2)≡[0,1]×[0,1]    dalla funzione

                                                                                         {({{ (Subscript[x, n+1]) = (Subscript[x, n]) + (Subscript[y, n]) + (ϵ/(2π)) sin(2π (Subscript[x, n]))    ( mod 1)}, { (Subscript[y, n+1])= (Subscript[y, n]) + (ϵ/(2π)) sin(2π (Subscript[x, n]))    ( mod 1)         }})

L'evoluzione di un punto iniziale ((Subscript[x, 0]),(Subscript[y, 0]))è data dalle successive applicazioni di tale funzione e l'insieme dei punti così ottenuti è detto orbita di ((Subscript[x, 0]),(Subscript[y, 0])).

Per ϵ=0 nei punti dell'orbita non cambia il valore iniziale di (Subscript[y, 0]) mentre l' ascissa aumenta ad ogni passo della stessa quantità.

Per ϵ>0 nei punti dell'orbita cambiano sia l'ascissa che l'ordinata ad ogni passo.

Facciamo degli esempi numerici

S M :  {({{(Subscript[x, n+1]) = (Subscript[x, n]) + (Subscript[y, n]) + (ϵ/(2π)) sin(2π (Subscript[x, n]))   (mod 1)}, {(Subscript[y, n+1]) = (Subscript[y, n]) + (ϵ/(2π)) sin(2π (Subscript[x, n]))             (mod 1)}})

f[eps_][{x_,y_}]:={Mod[x+y+eps/(2Pi)*Sin[2Pi*x],1],Mod[y+eps/(2Pi)*Sin[2Pi*x],1]}

orbita[x0_,y0_,eps_,nmax_]:=Block[{k,lista={},punto={x0,y0}},
For[k=1,k≤nmax,k++,
AppendTo[lista,punto];
punto=f[eps][punto];
];Return[lista];
]

orbita[0.1,0.618,0,4]

orbita[0.1,0.618,1,4]

ListPlot[orbita[0.1,0.618,0,4]]

ListPlot[orbita[0.1,0.618,0,4],PlotRange->{{0,1},{0,1}},AspectRatio->Automatic,PlotStyle->PointSize[0.05]]

Options[ListPlot]

disorbcorta[{x0_,y0_},eps_,nmax_]:=Block[{lista},lista=orbita[x0,y0,eps,nmax];
Return[ListPlot[lista,PlotRange->{{0,1},{0,1}},PlotStyle->PointSize[0.02],AxesLabel->{x,y},AspectRatio->Automatic,ImageSize->500,ColorFunction->Function[{x,y},Hue[x]],Epilog->{{PointSize[0.02],RGBColor[1,0,0],Point[{x0,y0}]},{PointSize[0.02],RGBColor[0,1,0],Point[lista[[-1]]]}}]];
]

Manipulate[disorbcorta[{0.1,.618},eps,n],{eps,0,1.5,.01,Appearance->"Open"},{n,1,100,1,Appearance->"Open"},ControlPlacement->Left]

Manipulate[disorbcorta[{x[[1]],x[[2]]},eps,n],{eps,0,1.5,.01,Appearance->"Open"},{n,1,100,1,Appearance->"Open"},{x,{0,0},{1,1},Appearance->"Open"},ControlPlacement->Left]

Manipulate[disorbcorta[{(Subscript[x, 0]),(Subscript[y, 0])},eps,n],{eps,0,1.5,.01,Appearance->"Open"},{n,1,100,1,Appearance->"Open"},{(Subscript[x, 0]),0,1,Appearance->"Open"},{(Subscript[y, 0]),0,1,Appearance->"Open"},ControlPlacement->Left]

S M :  {({{(Subscript[x, n+1]) = (Subscript[x, n]) + (Subscript[y, n]) + (ϵ/(2π)) sin(2π (Subscript[x, n]))   (mod 1)}, {(Subscript[y, n+1]) = (Subscript[y, n]) + (ϵ/(2π)) sin(2π (Subscript[x, n]))             (mod 1)}})

orbita[0.1,0.618,1,200]

disorblunga[{x0_,y0_},eps_,nmax_]:=Block[{lista},lista=orbita[x0,y0,eps,nmax];
Return[ListPlot[lista,PlotRange->{{0,1},{0,1}},PlotStyle->PointSize[0.001],AxesLabel->{x,y},AspectRatio->Automatic,ImageSize->500,Epilog->{{PointSize[0.02],RGBColor[1,0,0],Point[{x0,y0}]},{PointSize[0.02],RGBColor[0,1,0],Point[lista[[-1]]]}}]];
]

disorblunga[{x0_,y0_},eps_,nmax_]:=Block[{lista},lista=orbita[x0,y0,eps,nmax];
Return[ListPlot[lista,PlotRange->{{0,1},{0,1}},PlotStyle->PointSize[0.001],AxesLabel->{x,y},AspectRatio->Automatic,ImageSize->500,Epilog->{{PointSize[0.02],RGBColor[1,0,0],Point[{x0,y0}]},{PointSize[0.02],RGBColor[0,1,0],Point[lista[[-1]]]}}]];
]

disorblunga2[{x0_,y0_},eps_,nmax_,range_]:=Block[{lista},lista=orbita[x0,y0,eps,nmax];
Return[ListPlot[lista,PlotRange->range,PlotStyle->PointSize[0.001],AxesLabel->{x,y},AspectRatio->Automatic,ImageSize->500,Epilog->{{PointSize[0.02],RGBColor[1,0,0],Point[{x0,y0}]},{PointSize[0.02],RGBColor[0,1,0],Point[lista[[-1]]]}}]];
]

Manipulate[disorblunga[{0.1,.618},eps,n],{eps,0,1.5,.01,Appearance->"Open"},{n,1,10000,1,Appearance->"Open"},ControlPlacement->Left]

Manipulate[disorblunga[{x[[1]],x[[2]]},eps,n],{{eps,.7},0,1.5,.01,Appearance->"Open"},{{n,3000},1,10000,1,Appearance->"Open"},{{x,{.2,.7}},{0,0},{1,1},Appearance->"Open"},ControlPlacement->Left]

Manipulate[disorblunga2[{x[[1]],x[[2]]},eps,n,{{0,0.4},{0.5,0.6}}],{{eps,.7},0,1.5,.01,Appearance->"Open"},{{n,3000},1,10000,1,Appearance->"Open"},{{x,{.2,.7}},{0,0},{1,1},Appearance->"Open"},ControlPlacement->Left]

Manipulate[disorblunga[{(Subscript[x, 0]),(Subscript[y, 0])},eps,n],{{eps,.971},0,1.5,.01,Appearance->"Open"},{{n,5000},1,10000,1,Appearance->"Open"},{(Subscript[x, 0]),0,1,.01,Appearance->"Open"},{{(Subscript[y, 0]),.488},0,1,.01,Appearance->"Open"},ControlPlacement->Left]

cercaop[p_,q_,eps_]:=Block[{x=.5,y=.6,a=0,b=1,ls,xf=0,yf,ysum,err,res},
While[Abs[x-xf]>10^(-8),
ls=orbita[x,y,eps,q+1];
xf=ls[[q+1,1]];
yf=ls[[q+1,2]];
ls=Drop[ls,-1];
ysum=Apply[Plus,Map[#[[2]]&,ls]];
If[Abs[ysum-p]<10^(-8),Return[AppendTo[ls,{xf,yf}]]];
If[ysum>p,b=y,a=y];
y=(a+b)/2];
Return[AppendTo[ls,{xf,yf}]];
];

Manipulate[ListPlot[cercaop[p,q,eps]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.01],RGBColor[1,0,0]},AspectRatio->Automatic],{{p,5},1,21,1,Appearance->"Open"},{{q,8},1,21,1,Appearance->"Open"},{eps,.92,2,Appearance->"Open"},ControlPlacement->Left]

Manipulate[Show[ListPlot[Flatten[Join[Table[orbita[0.1,m/50,eps,1000],{m,1,49}],Table[orbita[0.5,m/50,eps,1000],{m,1,49}]],1]//Evaluate,PlotRange->{{0,1},{0,1}},ImageSize->500,AspectRatio->Automatic,AxesLabel->{x,y},PlotStyle->PointSize[0.001]],ListPlot[cercaop[p,q,eps]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.01],RGBColor[1,0,0]},AspectRatio->Automatic]],{{p,3},1,21,1,Appearance->"Open"},{{q,5},1,21,1,Appearance->"Open"},{{eps,.973},0,1.5,Appearance->"Open"},ControlPlacement->Left]

Show[ListPlot[Flatten[Join[Table[orbita[0.1,m/50,.973,1000],{m,15,40,.25}],Table[orbita[0.5,m/50,.973,1000],{m,15,40,.5}]],1]//Evaluate,PlotRange->{{0,1},{0,1}},ImageSize->700,AspectRatio->Automatic,AxesLabel->{x,y},PlotStyle->PointSize[0.001]],ListPlot[cercaop[1,2,.973]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.008],RGBColor[1,0,0]},AspectRatio->Automatic],ListPlot[cercaop[2,3,.973]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.008],RGBColor[1,0,0]},AspectRatio->Automatic],ListPlot[cercaop[3,5,.973]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.008],RGBColor[1,0,1]},AspectRatio->Automatic],ListPlot[cercaop[5,8,.973]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.008],RGBColor[0,1,0]},AspectRatio->Automatic],ListPlot[cercaop[8,13,.973]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.008],RGBColor[0,1,0]},AspectRatio->Automatic],ListPlot[cercaop[13,21,.973]//Evaluate,PlotRange->{{0,1},{0,1}},PlotStyle->{PointSize[0.008],RGBColor[1,0,0]},AspectRatio->Automatic]]