(* ______________________ birkhoff_HH.mth _________________________ In questo file sono raccolte le istruzioni che costituiscono il codice per effettuare (all'interno dell'ambiente di programmazione messo a disposizione dal software Mathematica) la costruzione della forma normale di Birkhoff per il modello di Henon-Heiles, nel caso in cui siano non-risonanti le velocita' angolari w1 e w2, corrispondenti al limite di piccole oscillazioni. *) (* Questo codice e' stato realizzato per la prima volta, in forma simile a quella presente, dalla Dr.ssa Chiara Caracciolo in occasione della preparazione della sua tesi di laurea, svolta sotto la supervisione del prof. Locatelli e terminata con la discussione finale nel corso della sessione di settembre 2016. Le successive variazioni sono state apportate con il permesso dell'autrice della versione originaria. Ultima modifica effettuata in data 02/05/2019. *) (* Definizioni iniziali -- INIZIO -- *) (* Il parametro R1 regola il troncamento degli sviluppi della Hamiltoniana rispetto al grado complessivo nelle variabili canoniche (x,y) o in quelle corrispondenti complesse (z, izb). *) R1 = 5; (* fs e' una lista, i cui elementi non sono altro che i termini polinomiali omogenei che compaiono nell'espansione delle Hamiltoniane. Ad esempio, al termine delle tre istruzioni seguenti, fs[[1]] contiene l'espansione del termine cubico di $H^{(0)}$ (in notazione TeX) espresso in funzione delle coordinate canoniche complesse (z, izb); fs[[2]] include l'espansione del termine quartico, etc. *) cubterm = x1^2 * x2 - x2^3 / 3; fs = Table[0, {l,1,R1-1}]; fs[[1]] = Expand[cubterm/.{x1->(z1-I*izb1)/N[Sqrt[2],16],x2->(z2-I*izb2)/N[Sqrt[2],16]}]; (* Definizione dei valori delle velocita' angolari w1 e w2. *) w1 = 1; w2 = N[(Sqrt[5] - 1) / 2, 16]; (* Definizione di Znf0, che include la parte iniziale gia' in forma normale, cioe' la coppia di oscillatori armonici. *) Znf0 = -I * w1 * z1 * izb1 -I * w2 * z2 * izb2; (* Definizione iniziale della lista chi, i cui elementi, alla fine del processo di normalizzazione, conterranno le espansioni delle funzioni generatrici $\chi^{(1)}$, $\chi^{(2)}$, etc. *) chi=Table[0, {l,1,R1-1}]; (* Definizioni iniziali -- FINE -- *) (* Il seguente module "poisson" restituisce l'espansione della parentesi di Poisson {f,g}, dove f e g sono le due funzioni in argomento, quando viene effettuata la chiamata di "poisson". *) poisson[f_, g_] := Module[{ris}, ris=0; ris = ris + D[f,z1] * D[g,izb1] - D[f,izb1] * D[g,z1]; ris = ris + D[f,z2] * D[g,izb2] - D[f,izb2] * D[g,z2]; ris]; (* Il seguente module "omosolv" restituisce l'espansione della generatrice X che risolve l'equazione omologica { Znf0 , X } + fr = Zr , dove fr e' il polinomio omogeneo che viene passato in argomento quando viene effettuata la chiamata di "omosolv"; inoltre, Zr altro non e' che "la media angolare" di fr, quando si adottano variabili di angolo azione. Equivalentemente, i monomi che compaiono nello sviluppo di Zr sono tali che l'esponente di z1 e' uguale a quello di izb1 e l'esponente di z2 e' uguale a quello di izb2. *) omosolv[fr_] := Module[{ris,frlist,nterms,j,pertterm,l1,l2,lb1,lb2,denom}, ris = 0; frlist = MonomialList[Expand[fr],{z1,z2,izb1,izb2}]; nterms = Length[frlist]; For[j=1,j<=nterms,j++, pertterm = frlist[[j]]; l1 = Exponent[pertterm, z1]; l2 = Exponent[pertterm, z2]; l1b = Exponent[pertterm, izb1]; l2b = Exponent[pertterm, izb2]; denom = I * ( (l1 - l1b) * w1 + (l2 - l2b) * w2 ); If[ (Abs[l1 - l1b] + Abs[l2 - l2b]) > 0, ris = ris - pertterm / denom; ]; ]; ris]; (* Il seguente module "prendicoefmax" non fa altro che restituire il coefficiente massimo (in valore assoluto) tra quelli che compaiono nello sviluppo dell'espressione expr, la quale e' funzione di (z1, z2, izb1, izb2). *) prendicoefmax[expr_] := Module[{ris,exprlist,nterms,j,monomio, coef,l1,l2,lb1,lb2}, ris = 0; exprlist = MonomialList[Expand[expr], {z1,z2,izb1,izb2}]; nterms = Length[exprlist]; For[j=1,j<=nterms,j++, monomio = exprlist[[j]]; If[Length[monomio] > 0, coef = monomio[[1]]; If[ Abs[coef] > ris, ris = Abs[coef]; ]; ]; ]; ris]; (* Il seguente module "filtra" restituisce un'espansione troncata rispetto a quella dell'espressione expr, la quale e' funzione di (z1, z2, izb1, izb2), in modo tale da omettere tutti i termini, i cui coefficienti (in valore assoluto) sono minori o uguali del valore soglia. *) filtra[expr_, soglia_] := Module[{ris,exprlist,nterms,j,monomio, coef,l1,l2,lb1,lb2}, ris = 0; exprlist = MonomialList[Expand[expr], {z1,z2,izb1,izb2}]; nterms = Length[exprlist]; For[j=1,j<=nterms,j++, monomio = exprlist[[j]]; If[Length[monomio] > 0, coef = monomio[[1]]; If[ Abs[coef] > soglia, ris = ris + monomio; ]; ]; ]; ris]; (* Il seguente module "passonorm" restituisce l'espansione della Hamiltoniana $H^{(r)}$, dal termine cubico in poi, dopo che viene effettuato l'r-esimo passo di normalizzazione. L'espansione della Hamiltoniana $H^{(r-1)}$ (cioe' relativa all'r-1-esimo passo di normalizzazione) e' rappresentata dal termine Znf0 e dalla lista fs, che viene passata in argomento, assieme al valore dell'indice r, quando viene effettuata la chiamata di "passonorm". *) passonorm[fs_,r_] := Module[{ris,s,j,jf,fix,lieterm,summand}, ris = fs; For[s=1,saz1*I, z1->1,izb2->az2*I, z2->1}; For[r=1,raz1*I, z1->1,izb2->az2*I, z2->1}; ]; (* Calcolo esplicito delle espansioni delle velocita' angolari in funzione delle sole azioni az1 e az2. Ovviamente, il seguente calcolo viene effettuato sostituendo a tutta la Hamiltoniana $H^{(R1-1)}$ la sua sola parte di forma normale. *) w1= D[Znftot, az1]; w2= D[Znftot, az2]; (* Il seguente module "scriviomega" si occupa di scrivere i termini dell'espansione w su di un file esterno individuato da stream; sia w che stream vengono passati in argomento, quando viene effettuata la chiamata di "scriviomega". Ovviamente, si intende che w sia una velocita' angolare, la cui espansione e' funzione delle sole azioni az1 e az2. *) scriviomega[w_, stream_] := Module[{wnew,wlist,numterm,monomio,i, esp1,esp2,coef,recoef,imcoef}, wnew = Expand[Simplify[w]]; wlist = MonomialList[wnew]; numterm = Length[wlist]; WriteString[stream, numterm, "\n"]; For[i=1, i<=numterm, i++, monomio= wlist[[i]]; esp1 = Exponent[monomio,az1]; esp2 = Exponent[monomio,az2]; coef = monomio/(az1^esp1*az2^esp2); recoef= N[Re[coef],16]; imcoef= N[Im[coef],16]; WriteString[stream,esp1,"\t",esp2,"\t", recoef,"\t","\t" imcoef,"\n"]; ]; WriteString[stream, "\n\n"]; ]; (* Apertura del file di output su cui successivamente vengono scritte le espansioni delle velocita' angolari w1 e w2, grazie a due oppotune chiamate del module "scriviomega". *) s=OpenWrite["omega.out"]; scriviomega[w1, s]; scriviomega[w2, s]; (* Chiusura del file di output *) Close[s];