Einführung in Mathematica XIV

ϵ-Algorthmus

© Jürgen Schnack, 2009/05/15

Problem

Habe eine Folge, deren Grenzwert ich aus wenigen Anfangsgliedern mit guter Genauigkeit extrapolieren möchte.

Beispiel: Berechne die Grundzustandsenergie für antiferromagnetische Spinringe für N=4, 6, 8, ..., 26. Extrapoliere zu N-> ∞.
Finite Size Scaling.

Echte Simulationsdaten
Spinring, af Kopplung, s=1/2

E1=Table[0.0,{m,1,26}];

E1[[2]]=-3.0;
E1[[3]]=-1.5;
E1[[4]]=-4.0;
E1[[5]]=-3.7360679775;
E1[[6]]=-5.6055512755;
E1[[7]]=-5.7103585137;
E1[[8]]=-7.3021868179;
E1[[9]]=-7.5945995679;
E1[[10]]=-9.0308927090;
E1[[11]]=-9.4378727250;
E1[[12]]=-10.7747818349;
E1[[13]]=-11.2591686595;
E1[[14]]=-12.5270990671;
E1[[15]]=-13.0673351449;
E1[[16]]=-14.2845927212;
E1[[17]]=-14.8670403306;
E1[[18]]=-16.0454981741;
E1[[19]]=-16.6609769085;
E1[[20]]=-17.8087730598;
E1[[21]]=-18.4508059644;
E1[[22]]=-19.5737613035;
E1[[24]]=-21.3400290294;
E1[[26]]=-23.107277686;

EProN=Table[{1/i,-E1[[i]]/i},{i,2,Length[E1]}];

ListPlot[EProN,PlotRange->{{0,0.55},{0.40,1.60}},Frame->True,PlotStyle->{PointSize[0.03],Red},BaseStyle->{FontFamily->Helvetica,FontSize->18},AxesLabel->{"1/i","E/i"},PlotLabel->"Energieeigenwerte"]

Merke

Manchmal ist nicht die ganze Folge zu gebrauchen, sondern nur eine Teilfolge.

Hier: Es gibt offensichtlich Unterschiede zwischen geradzahligen und ungeradzahligen Ringen (odd even effect). Der Grenzwert ist zwar gleich, aber für Beschleunigungsverfahren ist es besser, nur eine Teilfolge, z.B. die der geraden N zu nehmen.

Teilfolge

EG1=Table[-E1[[i]]/i,{i,2,Length[E1],2}]

ListPlot[EG1,Frame->True,PlotRange->{0.80,1.60},PlotStyle->{PointSize[0.03],Red},BaseStyle->{FontFamily->Helvetica,FontSize->18},AxesLabel->{"i","E/i"},PlotLabel->"Energieeigenwerte"]

Raten Sie den Grenzwert!
Stellen Sie sich vor, welchen Unfug eine Polynomapproximation ergeben würde!

ϵ-Feld anlegen (Indizes!!!)

kmax=Length[EG1]-1; (* von 0 bis 13 *)
pmax=kmax+1;
ϵ=Table[0.0,{p,-1,pmax},{k,0,kmax}];

ϵ[[p+2,k+1]] !

For[k=0,k≤kmax,k++,
ϵ[[0+2,k+1]]=EG1[[k+1]]];
Print[TableForm[ϵ]]

1. Generation: p=1

p=1;
For[k=0,k≤kmax-1,k++,
ϵ[[p+2,k+1]]=ϵ[[p-2+2,k+1+1]]+
1.0/(ϵ[[p-1+2,k+1+1]]-ϵ[[p-1+2,k+1]])];
Print[TableForm[ϵ]]

Alle weiteren Generationen

For[p=2,p≤pmax-1,p++,
For[k=0,k≤kmax-p,k++,
ϵ[[p+2,k+1]]=ϵ[[p-2+2,k+1+1]]+
1.0/(ϵ[[p-1+2,k+1+1]]-ϵ[[p-1+2,k+1]])];
Print[TableForm[ϵ]];
Print["----------------------------------------------"]]

Zusammenfassung

Results=Table[ϵ[[p+2,0+1]],{p,0,pmax,2}]

Das wahre Ergebnis ist in diesem Fall (af Ring, s=1/2) bekannt und geht auf Hans Bethe (Bethe-Ansatz) zurück: 2 ln(2)-(1/2)=0.886294

N[2*Log[2.0]-1/2]

Die Extrapolation ist deutlich näher am wahren Ergebnis als das höchste Folgenglied (0.888741).

Levin-u- und Levin-t-Methode

EG1D=Table[EG1[[i]]-EG1[[i-1]],{i,1,Length[EG1]}];
EG1D[[1]]=EG1[[1]];

T[n_]:=Sum[(-1)^(m-1) m^(n-1) Binomial[n,m] EG1[[m]]/EG1D[[m]],{m,1,n}]/Sum[(-1)^(m-1) m^(n-1) Binomial[n,m] 1/EG1D[[m]],{m,1,n}];
U[n_]:=Sum[(-1)^(m-1) m^(n-2) Binomial[n,m] EG1[[m]]/EG1D[[m]],{m,1,n}]/Sum[(-1)^(m-1) m^(n-2) Binomial[n,m] 1/EG1D[[m]],{m,1,n}];

TList=Table[N[T[n],10],{n,1,Length[EG1]}]
UList=Table[N[U[n],10],{n,1,Length[EG1]}]

ϵ-Methode einfach und schnell.

Levin-Methoden ausgefeilter; konvergieren schneller und genauer. 
Vgl. wahres Ergebnis 2 ln(2)-(1/2)=0.886294