Oscillatore anarmonico
H= (1/2)(p^2) + (1/2)(x^2) + g(1/2)(x^4)

Parametri del problema

g = 10;     (* costante di anarmonicità *)

livelli = {1,2,10};

Primo Metodo

Si scrivono nella base dei polinomi di Hermite le matrici per
gli operatori  (x^2) , (x^4)  e   (p^2). Si costruisce l'Hamiltoniana, (1/2)p2 + (1/2)x2 + g(1/2)x4,
 e si diagonalizza la matrice.

Calcolo  degli operatori nella base dell'oscillatore  armonico

Conviene scrivere una funzione che restituisca le tre matrici,
in modo da poterla riutilizzare. Gli elementi di matrice si trovano nel testo.

matrici[dim_]:=
Module[{x2array,p2array,x4array,x1array,x2,p2,x4}, 
n = i-1;
x2array=Normal[SparseArray[{{i_,i_}->  (1/4)(2n+1),{i_,k_}/;k-i==2->  (1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}]];
x2 =(x2array+Transpose[x2array]);
p2array =Normal[SparseArray[{{i_,i_}->  (1/4)(2n+1),{i_,k_}/;k-i==2->  -(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}]];
p2 = (p2array+Transpose[p2array]);
x4array=Normal[SparseArray[{{i_,i_}->  (3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==2-> ((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==4-> (1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}]];
x4 =  (x4array+Transpose[x4array]);
{p2,x2,x4}
];

Prova:

{p2,x2,x4} = matrici[40];

ham =N[ (1/2)p2 + (1/2)x2 + g(1/2)x4];

Sort[Eigenvalues[ham]][[livelli]]

{1.224593426404839`,4.2995158786769005`,68.5628166735637`}

{p2,x2,x4} = matrici[80];ham =N[ (1/2)p2 + (1/2)x2 + g(1/2)x4];ham =N[ (1/2)p2 + (1/2)x2 + g(1/2)x4];Sort[Eigenvalues[ham]][[livelli]]

1.2245870362635476`,4.299501734368853`,48.08842968530925`}

Controllo della convergenza della procedura, aumentando il numero dei vettori di base. Per la stampa,prima si definisce la stringa vuota eigStampa, che viene via via riempita con la dimensione considerata e i 3 autovalori che si vogliono stampare.

eigStampa = {};
Do[{dim = j*20;
{p2,x2,x4} = matrici[dim];
ham =N[ (1/2)p2 + (1/2)x2 + g(1/2)x4];
eigH = Sort[Eigenvalues[ham]][[livelli]];
eigStampa = Append[eigStampa,Prepend[eigH,dim]];
},{j,1,7}]

TableForm è un modo comodo di stampare una tabella.
TableHeadings  definisce la prima riga.

NumberForm[TableForm[eigStampa,TableHeadings-> {None,{"base","  E0 ", "  E1 ","  E2" }}],{14,10}]

Valore abbstanza preciso per g=10  E = 1.2245870361 (stato fond.)

Secondo Metodo

Gli operatori hanno elementi di matrice solo fra stati con la stessa parità quindi conviene costruire matrici solo per
gli stati pari o solo per quelli dispari, come nella funzione seguente

Remove[p2,x2,x4];

Matrici in base dell'oscillatore armonico

matricipar[dimloc_,parita_]:=
Module[{x2array,p2array,x4array,x2,p2,x4,q2,y2,y4}, 
Clear[n,i,j];
n =2i-2+parita;
x2array=Normal[SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,j_}/;j-i==1->(1/2)(Sqrt[(n+1)(n+2)])},{dimloc,dimloc}]];
x2 = N[(x2array+Transpose[x2array])];
p2array =Normal[SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,j_}/;j-i==1->-(1/2)(Sqrt[(n+1)(n+2)])},{dimloc,dimloc}]];
p2 = N[(p2array+Transpose[p2array])];
x4array=Normal[SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,j_}/;j-i==1->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,j_}/;j-i==2->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dimloc,dimloc}]];
x4 =  N[(x4array+Transpose[x4array])];
{p2,x2,x4}
]

Il seguito è identico a prima: per esempio, per gli stati pari  (per i dispari, scrivere:  matricipar[dim,1] )

eigStampa = {};
Do[{dim = j*20;
{p2,x2,x4} = matricipar[dim,0];
ham =N[ (1/2)p2 + (1/2)x2 + g(1/2)x4];
eigH = Sort[Eigenvalues[ham]][[livelli]];
eigStampa = Append[eigStampa,Prepend[eigH,dim]];
},{j,1,5}]

NumberForm[TableForm[eigStampa,TableHeadings-> {None,{"base","  E0 ", "  E2 ","  E4" }}],{14,10}]

Si nota che la convergenza è molto più veloce che nel metodo 1.

{p2,x2,x4}=matricipar[20,0];
ham2 =N[ (1/2)p2 + (1/2)x2 + g(1/2)x4];
Sort[Eigenvalues[ham2]][[livelli]]

MatrixForm[ham2[[livelli,livelli]]]