High order expansion p.t.

Iterative algorithm (Sakurai-Onofri-Destri)

SetDirectory["c:/home/enrico/work/TBA/Math"];

mu = 4;
$IterationLimit = 4096*4

Code by R. De Pietri

SetAttributes[{Operatore}, HoldFirst];

Operatore/: Operatore[op_][x_+y_]:=Operatore[op][x]+Operatore[op][y]
Operatore/: Operatore[op_][ x_ + y_ ] := Operatore[op][x] + Operatore[op][y]
Operatore/: Operatore[op_][x_]  := 0                        /; FreeQ[x,Vettore] 
Operatore/: Operatore[op_][ x_  y_  ] := x Operatore[op][y] /; FreeQ[x,Vettore] 
Operatore/: x_ Operatore[op_][Vettore[y_,z_]] := 
               Operatore[op][Vettore[x y,z]] /; FreeQ[x,Vettore]

Vettore/: Vettore[c1_,x_] + Vettore[c2_,x_] := Vettore[ c1 + c2 , x]
Vettore/: c1_ Vettore[c2_,x_]               := Vettore[ c1  c2  , x]

Format[Vettore[x_,y_]] := SequenceForm[x,"|",y,">"]

sc[ x_ , y1_ + y2_ Vettore[ z_]] := sc[x,y1] + y2 sc[x, Vettore[z]]
sc[ y1_ + y2_ Vettore[ z_] , x ] := sc[y1,x] + y2 sc[Vettore[z], x]

sc[ x_ + y_, z_] := sc[x,z] + sc[y,z]
sc[ z_,x_ +  y_] := sc[z,x] + sc[z,y]


sc[ x_, y_]  := 0 /; FreeQ[y,Vettore] 
sc[ y_, x_]  := 0 /; FreeQ[y,Vettore] 

sc[ x_ y_ ,  z_]  := x sc[y,z] /; FreeQ[x,Vettore] 
sc[ z_, x_   y_]  := x sc[z,y] /; FreeQ[x,Vettore]

sc[ Vettore[x_,n_] , Vettore[y_,n_] ] := (x y/ (2 n + mu + 1))
sc[ Vettore[_ ,n_] , Vettore[_ ,m_] ] := 0 /; 
                     IntegerQ[ Simplify[ n - m ] ] && ( Simplify[n-m] != 0 )

a  = Operatore[a]  ;   
ad = Operatore[ad]  ;    
V  = Operatore[V]    ;   
Q   = Operatore[Q];

Operatore/: a[Vettore[x_,0]]     := 0
Operatore/: a[x_] := 0  /; FreeQ[x,Vettore] 
Operatore/: a[Vettore[c1_,n_]]    := Vettore[2 n^2 (mu + n)^2/((mu+2n)(2n+mu+1)) c1, n-1]

Operatore/: ad[x_] := 0 /; FreeQ[x,Vettore] 
Operatore/: ad[Vettore[c1_,n_]]    := Vettore[2(n+1)^2 (n+mu+1)^2 / ((2n+mu+1)(2n+mu+2)) *c1, n+1]

Operatore/: Q[x_] := 0 /; FreeQ[x,Vettore] 
Operatore/: Q[Vettore[c1_,n_]]    := Vettore[mu(mu(1+mu)+(1+mu)(2+mu)n+(2+mu)n^2)/(2n+mu)/(2n+mu+2) *c1, n]

Operatore[R[n_]][Vettore[_,n_]] = 0 ;
Operatore[R[m_]][Vettore[c1_,n_]] := Vettore[c1 1/(autovH0[n]-autovH0[m]) ,n] ;

Inizializza[ iteratore ] := Block[{}
,ClearAll[ iterazione , phi , delta ]
;phi[m_][0] := Vettore[1,m](2m + mu + 1)
;iterazione[m_,0] := True
    ;iterazione[m_,n_]:= Block[ {},
        iterazione[m,n-1];
        delta[m,n] =Simplify[ sc[ Vettore[1,m] , V @ phi[m][n-1] ] ] ;
      
        phi[m][n] = Simplify[Expand[
                    - Operatore[R[m]] @ V @ phi[m][n-1]
                    + Sum[ delta[m,k] Operatore[R[m]] @ phi[m][n-k] ,{k,1,n-1}] 
        ]] ;
        iterazione[m,n] = True ;
        Print["got correction ",n," state |",m,">"];
    ]
]

FUNZIONI DI CALCOLO DISPONIBILI

autovalore[m_,n_] := Block[ {} , 
                     iterazione[m,n] ; 
                     autovH0[m] + Sum[ delta[m,k] , {k,1,n}]  ]

autovettore[m_,n_] := Block[ {} , 
                      iterazione[m,n] ; 
                      Sum[ phi[m][k] , {k,1,n}]  ]

Principio di funzionamento del programma
  sia  H =  H_0 + p V  con  H_0 |n> = autovH0[n] |n>

  cerchiamo una soluzione di H |aut(n)> = Vaut(n) |aut(n)>
  in serie di potenze in P
  |aut(n)> = |n> + p |n,1> + p^2 |n,2> + p^3 |n,3> + ..
    Vaut(n)  = autovH0[n] + p delta[n,1] + p^2 delta[n,2] + ...

  allora si ha che la soluzione del problema e` data dalle  regole ricorsive

  delta[m,k] = <m| V |m,k-1>

  |m,k>  = R[m] ( V |m,k-1> + delta[m,1]   |m,k-1> + ... 
                            + delta[m,k-1] |m,1> )

Problema[1] := Block[{}

   ,Inizializza[iteratore]
   ;autovH0[ n_ ] := n(n+mu+1)
   ;V  @ x_ := eps * Q @ x - eps * ad @ x - eps *  a @ x

]

Problema[1]

ClearAll[q]

q[n_, 
          ord_] := q[n, 
            ord] = (1 + mu)/2  +  autovalore[n, ord]/eps/
                2 /. {eps -> x};

κ[n_, ord_] := 
        Series[(24  q[n, ord] - 6 (1 + 2 mu) (w + 1))/(1 - 
                  w) /. w -> 1/2 (1 - 1/x), {x, 0, 
            ord}];

This is the 2^nd order PT done by hand:

Eig[j_, m_] := 
      j (j + 1) - m/2 (m/2 + 1) (1 - x) - 
        x m^4/(16  j (j + 1)) + 
        x^2 ((4  j^2 - m^2)^4/(512 
                     j^3 (4  j^2 - 1)) - (4 (j + 1)^2 - 
                      m^2)^4/(512 (j + 1)^3 (4 (j + 1)^2 - 
                        1)))

Eig0[
          m_] := (j (j + 1) - m/2 (m/2 + 1) (1 - x) - 
              x m^4/(16  j (j + 1)) + 
              x^2 (-(4 (j + 1)^2 - 
                            m^2)^4/(512 (j + 
                              1)^3 (4 (j + 1)^2 - 1)))) /. 
          j -> m/2;

Myq[j_, m_] := (1 + m)/2 + Eig[j, m]/x/2;

Myκ[j_, m_] := 
        Series[(24  Myq[j, m] - 6 (1 + 2 m) (w + 1))/(1 - 
                  w) /. w -> 1/2 (1 - 1/x), {x, 0, 
            2}];

K[j_, m_] := 
        Simplify[
          Series[(24  Eig[j, m] + 6 (1 + 2  m) + 
                  6  x (1 - 2  m))/(1 + x), {x, 0, 
              2}]];

K2[j_, m_] :=  
        6 (1 + 4  j^2 + 4  j - m^2) +  
          x (12  m^2 - 24  j (j + 1) - 
                3/2  m^4/(j (j + 1))) +  
          x^2 (24 (j (j + 1) - m^2/2) + 
                3/2  m^4/(j (j + 1)) + 
                3 ((j^2 - 
                              m^2/4)^4/(j^3 (j^2 - 
                                1/4)) - ((j + 1)^2 - 
                              m^2/4)^4/((j + 1)^3 ((j + 1)^2 - 
                                1/4))));

K20[j_, m_] :=  
        6 (1 + 4  j^2 + 4  j - m^2) +  
          x (12  m^2 - 24  j (j + 1) - 
                3/2  m^4/(j (j + 1))) +  
          x^2 (24 (j (j + 1) - m^2/2) + 
                3/2  m^4/(j (j + 1)) + 
                3 (-((j + 1)^2 - 
                              m^2/4)^4/((j + 
                                1)^3 ((j + 1)^2 - 
                              1/4))));

K[j, 0]

K2[j, 0] - K[j, 0]

Table[K[j, m] - K2[j, m], {j, 2, 10}, {m, 1, 2}]

K[j, m] - K2[j, m]

Series[K2[j, m] /. x -> y/(2 - y), {y, 0, 2}]

κ[0, 10] /. {x -> λ/(2 - λ)}

Series[Normal[%], {λ, 0, 10}]

% // N

k03 = CoefficientList[%, λ];

DumpSave["kappa10"];

Sum[k03[[j + 1]] (λ^j), {j, 0, 10}] /. λ -> 
        1

N[%, 20]

6 ((mu^2) - (((mu - 1))^2))

42.00100251920617356075019678/42 - 1

ListPlot[
      Table[Log[(Sum[
                  k03[[j + 1]] (λ^j), {j, 0, 
                    n}] /. λ -> 1) - 42], {n, 1, 100}], 
      PlotJoined -> True]