Bacini multipli di attrazione: il metodo di Newton

Un metodo molto famoso per trovare zeri di funzioni è il metodo di Newton. Questo metodo si basa sulla convergenza veloce delle successione ... agli zeri della funzione.

Questo metodo è un algoritmo di ricerca degli zeri di una funzione (complessa) f[z], e si basa sull'approssimazioe lineare della funzione

Supponiamo data una funzione f[z] ed ipotizziamo che lo zero che ci interessa si trovi in (Subscript[z, 0]).

Sappiamo che vale f[z]=f[(Subscript[z, 0])]+f'[(Subscript[z, 0])] (z-(Subscript[z, 0])) + O((((z-(Subscript[z, o])))^2)). Approssimiamo f[z] con (Subscript[g, (Subscript[z, 0])])[z]=f[(Subscript[z, 0])]+f'[(Subscript[z, 0])] (z-(Subscript[z, 0])).

Lo zero di (Subscript[g, (Subscript[z, 0])])[z] sarà (Subscript[z, 1])=(Subscript[z, 0])-((f[(Subscript[z, 0])])/(f'[(Subscript[z, 0])]))

Si ottiene così una successione di punti (Subscript[z, 0]),(Subscript[z, 1]),... che dovrebbe convergere verso uno zero di f[z]

Potremmo anche mettere in piedi un metodo che si basa sull'approssimazioe al secondo ordine della funzione f[z]

Supponiamo data una funzione f[z] ed ipotizziamo che lo zero che ci interessa si trovi in (Subscript[z, 0]).

Sappiamo che vale f[z]=f[(Subscript[z, 0])]+f'[(Subscript[z, 0])] (z-(Subscript[z, 0])) +((f''[(Subscript[z, 0])])/2)(((z-(Subscript[z, o])))^2)+O((((z-(Subscript[z, o])))^3)). Approssimiamo f[z] con (Subscript[g, (Subscript[z, 0])])[z]=f[(Subscript[z, 0])]+f'[(Subscript[z, 0])] (z-(Subscript[z, 0]))+((f''[(Subscript[z, 0])])/2)(((z-(Subscript[z, o])))^2).

z/.Solve[f+fp(z-w)+fpp/2(z-w)^2==0,z]

Si ottengono quindi 2 soluzioni (Subscript[z, 11]),(Subscript[z, 12]), che sono -((f'[(Subscript[z, 0])]±(Sqrt[((f'[(Subscript[z, 0])])^2)-2 f[(Subscript[z, 0])]f''[(Subscript[z, 0])]])+f''[(Subscript[z, 0])](Subscript[z, 0]))/(f''[(Subscript[z, 0])]))

Osserviamo che, se effettivamente f[(Subscript[z, 0])]=0, allora i due punti ottenuti sono

Simplify[{((f'[(Subscript[z, 0])]+(Sqrt[((f'[(Subscript[z, 0])])^2)-2 f[(Subscript[z, 0])] f''[(Subscript[z, 0])]])+f''[(Subscript[z, 0])] (Subscript[z, 0]))/(f''[(Subscript[z, 0])])),((f'[(Subscript[z, 0])]-(Sqrt[((f'[(Subscript[z, 0])])^2)-2 f[(Subscript[z, 0])] f''[(Subscript[z, 0])]])+f''[(Subscript[z, 0])] (Subscript[z, 0]))/(f''[(Subscript[z, 0])]))}/.f[(Subscript[z, 0])]->0,f'[(Subscript[z, 0])]<0]

sdoppia[w_]:=N[{((f'[w]+(Sqrt[((f'[w])^2)-2 f[w]f''[w]])+f''[w]w)/(f''[w])),((f'[w]-(Sqrt[((f'[w])^2)-2 f[w]f''[w]])+f''[w]w)/(f''[w]))}]
vai[lista_]:=Flatten[Map[sdoppia,lista]]

z/.NSolve[f[z]==0,z]

f[z_]:=(z^3)-1

sdoppia[1]

approx[n_]:=Nest[vai,{0.1},n]

approx[5]

Non sembra convergere a nulla

Esercizio: disegnare i bacini di attrazione della iterazione di Newton per i polinomi ciclotomici (z^n)-1

Secondo tentativo, funziona meglio per i polinomi ciclotomici

cambia[x_]:=If[Head[x]===Real,x,1]

colore[z0_,F_]:=cambia[Arg[Nest[F,N[z0],50]]/2π]

Con le radici terze

Gli zeri veri della funzione sono

zeri=Map[Extract[#,{1,2}]&,N[Solve[(z^3)-1==0,z]]];

La funzione da iterare è (Subscript[z, 0])-((f[(Subscript[z, 0])])/(f'[(Subscript[z, 0])]))

z-(((z^3)-1)/(D[(z^3)-1,z]))//Simplify

F[z_]:=(1/(3 (z^2)))+((2 z)/3)

Creiamo una tabella in cui ogni entrata

tabella=Table[colore[a+b I,F],{b,-1,1,0.01},{a,-1,1,0.01}];

Flatten[tabella]//Union

ListContourPlot[tabella]

con le radici quinte, abbastanza dettagliato

zeri=Map[Extract[#,{1,2}]&,N[Solve[(z^5)-1==0,z]]];
F[z_]:=(1/(5 (z^4)))+((4 z)/5)

tabella=Table[colore[a+b I,F],{b,-2,2,0.01},{a,-2,2,0.01}];

ListContourPlot[tabella]

z-((z^8)-1)/D[z^8-1,z]//Simplify

Con la funzione seno

zeri={0,-π,π}
F[z_]:=z-((Sin[z])/(Cos[z]))

tabella=Table[colore[a+b I,F],{b,-0.2,0.2,0.001},{a,(π/3),(2/3)π,0.001}];

ListContourPlot[tabella]

Laghi di Wada

Consideriamo il sistema dinamico del pendolo forzato–smorzato.

Un pendolo di lunghezza ℓ e massa m, soggetto alla forza di gravità, soddisfa all'equazione differenziale

m ℓ (x^(••))=-γ (x^•)- m g Sin[x]

Aggiungere un termine forzante vuol dire aggiungere all'equazione un termine

m ℓ (x^(••))=-γ (x^•)- m g Sin[x]+f[t]

dove f[t]:ℝ->ℝ è detto termine forzante (o semplicemente forzante).

Cambiando l'unità di misura dello spazio, del tempo e della massa, si a–dimensionalizza l'equazione e ci si riduce a

(x^(••))=-γ (x^•)-Sin[x]+f[t]

Il campo vettoriale associato a questa equazione differenziale è non–autonomo (ovvero dipende da t)

X=(({{v}, {-γ v - Sin[x]+f[t]}}))

L'unicità delle soluzioni, in questo caso, è vera prendendo in considerazione il campo vettoriale in (Subscript[ℝ, x])×(Subscript[ℝ, v])×(Subscript[ℝ, t]) le cui componenti sono

X=(({{v}, {-γ v - Sin[x]+f[t]}, {1}}))

Ovviamente, un tale campo vettoriale non può avere soluzioni di equilibrio (non è mai nullo). Ma potrebbe avere soluzioni periodiche.

In particolare, quando il termine forzante è periodico ad esempio f[t]=Cos[t], allora

Prop. le soluzioni periodiche hanno periodo che è multiplo del periodo di f

Esercizio P1. risolvere il sistema e mostrare che le soluzioni tendono verso soluzioni periodiche. Rappresentare tali soluzioni periodiche attraverso grafici ed animazioni

Sembra insomma che una buona parte dei dati iniziali convergano verso un'orbita periodica

Volendo determinare verso quale orbita periodica venga attratto un dato punto iniziale, dovremmo

0. scegliere un'orbita periodica, e determinare una palla all'interno della quale tutti i dati iniziali convergono verso tale orbita. Una volta scelta tale palla, le altre palle sono le 2k π traslate di tale palla 
2. dato un punto qualsiasi, calcolare la sua orbita fino a che questa non sia caduto all'interno di una di tali palle. Assegnare a quel punto il valore k
3. disegnare i punti del piano in accordo con il valore di k assegnato ad essi

Svolgiamo l' esercizio usando l'integratore di Mathematica

Questo calcola la soluzione

sol[{x0_,v0_}]:=Module[{},
{x[t],v[t]}/.NDSolve[{x'[t]==v[t],v'[t]==Cos[t]-0.1 v[t]-Sin[x[t]],x[0]==x0,v[0]==v0},{x[t],v[t]},{t,0,250}][[1]]]

Questo produce l' andamento della x come funzione di t

disegna[{x0_,v0_}]:=Module[{},
s=sol[{x0,v0}][[1]];
Plot[s,{t,0,250},PlotRange->{9π,-9π},Ticks->{Automatic,Table[π i, {i,-10,10}]},PlotStyle->{Hue[3v0],Thickness[((v0-1.7)/100)]}]]

Questo produce il filmato della soluzione

filmsom[{x0_,v0_}]:=Module[{},
s=sol[{x0,v0}];
fotogrammi[xmin_]:=ParametricPlot[s,{t,xmin,xmin+1},PlotRange->{{-5π,5π},{-2π,2π}},Ticks->{Table[π i, {i,-10,10}],Automatic}];
Animate[fotogrammi[xmin],{xmin,0,248,0.01},AnimationRate->2]]

Disegnamo l' andamento delle x, e vediamo che tendono a stabilizzarsi (dopo un po' di giri) ad un moto periodico

Show[Table[disegna[{0,i}],{i,1.8,2.1,0.05}]]

Visualizziamo tali moti periodici

filmsom[{0,1.8}]

Tendono a stabilizzarsi ... non tutti [13]

Scriviamo la mappa di Poincaré

Φ[{x_,v_}]:=sol[{x,v}]/.t->2π
Φ[n_][{x_,v_}]:=NestList[Φ,{x,v},n]

Vediamo che ogni dato iniziale vicino all'origine si stabilizza su di un particolare punto di un' orbita periodica, che è quella che dice Hubbard a pag 745

Last[Φ[100][{0,0}]]

Mostriamo come molti dati iniziali convergono verso la stessa soluzione

Show[
ListPlot[Φ[100][{0,0.02}],PlotStyle->Red],
ListPlot[Φ[100][{0,0.03}],PlotStyle->Green],
ListPlot[Φ[100][{0,0.04}],PlotStyle->Blue],
ListPlot[Φ[100][{0,0.01}],PlotStyle->Yellow],
ListPlot[Φ[100][{0,0}],PlotStyle->Magenta]]

Se sposto il dato iniziale di 2 π, ottengo l'orbita periodica spostata di 2π

Last[Φ[100][{2π,0.02}]]-2π

Vediamo il filmato dell' orbita periodica "base"

filmsom[{-2.046291629029782`,0.39269675159911344`}]

E di quella traslata di 2 π

filmsom[{-2.046291629029782`+ 2π,0.39269675159911344`}]

Vanno tutti verso lo stesso equilibrio

Coloriamo ora il disegno

(4/768)//N

solT[{x0_,v0_}]:=Module[{},
{x[t],v[t]}/.NDSolve[{x'[t]==v[t],v'[t]==Cos[t]-0.1 v[t]-Sin[x[t]],x[0]==x0,v[0]==v0},{x[t],v[t]},{t,0,60π}][[1]]]/.t->60π

array=Table[{i,j},{j,-2.,2.,0.0052083333333334},{i,-N[3π],N[2π],0.0153398078788565}];

arrayT=Map[solT,array,{2}];

Associamo ai punti raggiunti un valore di colore

col[{x_,v_}]:=Which[
(((x+2.046291506120066-2π))^2)+(((v-0.39269675159911344))^2)<0.1,0,
(((x+2.046291506120066))^2)+(((v-0.39269675159911344))^2)<0.1,1,
(((x+2.046291506120066+2π))^2)+(((v-0.39269675159911344))^2)<0.1,2,
1==1,3]

ArrayPlot[Map[col,arrayT,{2}],
ColorRules->{0->Red,1->Green,2->Blue,3->White}]

Dimensions[arrayT]

col[arrayT[[4,5]]]

Scriviamo un integratore

Proviamo con un Strang... studiamo l'appartenzenza in un modo brutale

Dobbiamo necessariamente scrivere un nostro integratore. Il campo vettoriale è

X=(({{v}, {-γ v - Sin[x]+f[t]}, {1}}))=(({{v}, {-γ v}, {1}}))+(({{0}, {- Sin[x]+f[t]}, {0}}))

I flussi dei due campi sono (Subscript[Φ, 1])[h][x,v,t]=(({{x+(v/γ)(1-(e^(-γ h)))}, {v (e^(-γ h))}, {t+h}})), (Subscript[Φ, 2])[h][x,v,t]=(({{x}, {v+(f[t]-Sin[x])h}, {t}}))

Quindi, usando un metodo splitting, definiamo le funzioni F1, che fa avanzare il tempo di h ed F2, che non muove il tempo

F1[h_][{x_,v_,t_}]:={x+v(1-(E^(-h))),v (E^(-h)),t+h}
F2[h_][{x_,v_,t_}]:={x,v +(Cos[t]-Sin[x])h,t}

Componiamole al tempo h

F2[h][F1[h][x,v,t]]//Simplify

Si ottiene una funzione fa avanzare t di h

FLTcomp=Compile[{{z,_Real,1}},
{z[[2]]-(E^(-π/1000))z[[2]]+z[[1]],(E^(-π/1000)) z[[2]]+((π Cos[(π/1000)+z[[3]]])/1000)-((π Sin[z[[2]]-(E^(-π/1000)) z[[2]]+z[[1]]])/1000),(π/1000)+z[[3]]}];

Prendiamo h=(π/100), così che si devono fare 200 iterazioni per arrivare a 2π, e quindi 9999 iterazioni per arrivare a 50π

FS[{x_,v_,t_}]:=Compile[{{z,_Real,1}},
F2[(1/2)(π/1000)][Nest[FLTcomp,F1[(π/1000)][F2[0.5(π/1000)][z]],99999]]]

array=Table[{i,j,0},{i,-N[3 π],N[3 π],0.1π },{j,-3.,3.,0.1}];

arrayFS=Map[FS,array,{2}]//Timing

col[{x_,v_,t_}]:=Which[
Abs[x+2.046291506120066-2π]<0.01,0,
Abs[x+2.046291506120066]<0.01,1,
Abs[x+2.046291506120066+2π]<0.01,2,
Abs[x+2.046291506120066+2π]≥0.01,3]

array[[2]]

%98[[2,40]]

Map[col,%98[[2]],{2}]

ArrayPlot[array,
ColorRules->{0->Red,1->Green,2->Blue,3->White}]

Proviamo con NDSolve per tempi abbastanza lunghi

sol[{x0_,v0_}]:=Module[{},
soluz={x[t],v[t]}/.NDSolve[{x'[t]==v[t],v'[t]==Cos[t]-0.1 v[t]-Sin[x[t]],x[0]==x0,v[0]==v0},{x[t],v[t]},{t,0,60π}][[1]];
soluz/.t->60π]

array=Table[sol[{i,j}],{i,-N[4 π],N[4 π],0.02π },{j,-3.,3.,0.01}]

col[{x_,v_}]:=Which[
(*Se la x è vicina a quella del "punto fisso" del ciclo base -2π, si ha 1*)
Abs[x+2.046291506120066-2π]<0.01,0,
(*Se la x è vicina a quella del "punto fisso" del ciclo base, si ha 1*)
Abs[x+2.046291506120066]<0.01,1,
(*Se la x è vicina a quella del "punto fisso" del ciclo base+2π, si ha 1*)
Abs[x+2.046291506120066+2π]<0.01,2,
(*Quest'ultima linea restituisce 3 se non ha passato nessuno dei test precedenti*)
1==1,3]

ArrayPlot[Map[col,array,{2}],
ColorRules->{0->Red,1->Green,2->Blue,3->White}]

Proviamo con un while

Facciamo il disegno

F1[h_][{x_,v_,t_}]:={x+v(1-(E^(-h))),v (E^(-h)),t+h}
F2[h_][{x_,v_,t_}]:={x,v +(Cos[t]-Sin[x])h,t}

Componiamole al tempo h

F2[h][F1[h][{x,v,t}]]//Simplify

Si ottiene una funzione fa avanzare t di h

FLTcomp=Compile[{{z,_Real,1}},
{z[[2]]-(E^(-π/1000))z[[2]]+z[[1]],(E^(-π/1000)) z[[2]]+((π Cos[(π/1000)+z[[3]]])/1000)-((π Sin[z[[2]]-(E^(-π/1000)) z[[2]]+z[[1]]])/1000),(π/1000)+z[[3]]}];

Prendiamo h=(π/100), così che si devono fare 200 iterazioni per arrivare a 2π, e quindi 9999 iterazioni per arrivare a 50π

F2[(1/2)(π/1000)][Nest[FLTcomp,F1[(π/1000)][F2[0.5(π/1000)][{0,0,0}]],99999]]//Timing

Val[{x_,v_}]:=Which[
(((x0-x))^2)+(((v0-v))^2)<0.1,0,
(((x0+2 π-x))^2)+(((v0-v))^2)<0.1,1,
(((x0-2 π-x))^2)+(((v0-v))^2)<0.1,-1,
1==1, False]

&&((((x0+2 π-a))^2)+(((v0-b))^2)>0.1)&&((((x0-2 π-a))^2)+(((v0-b))^2)>0.1)

FluisciFino[{x_,v_}]:=Module[{x0=-2.046291506120066,v0=0.39269675159911344,a=x,b=v,c=0},
While[((((x0-a-2π))^2)+(((v0-b))^2)>0.1),
Print[(((x0-a))^2)+(((v0-b))^2)]; {a,b,c}=FLTcomp[{a,b,c}]];
a]

FluisciFino[{0,0}]