(*Napisane przez Pawła Czachorowskiego. Wykorzystanie dowolne.*)
(*Wszelkie uwagi mile widziane.*)

(*Wariacje (w sensie potocznym, nie naukowym) na temat modelu dwóch stanów (Literatura: Dodatek D z "Idei chemii kwantowej L. Pieli" (Wydanie 2. Polskie 2011)).*)

(*Mamy układ, który może być opisany jako kombinacja liniowa dwóch stanów. Zakładamy ich normalizację. Właściwie matematycznie to taki sam problem, co w zadaniu z metodą Ritza i dwoma funkcjami Gaussa. Tyle, że tutaj nie zakładamy żadnej postaci funkcji (poza tym, że jest znormalizowana) i poszukujemy ogólnego rozwiązania oraz jego ciekawych przypadków.*)

(*Czynimy pewne założenia: funkcje (a zatem i elementy macierzowe H i S) są rzeczywiste, a także mówimy o stanach związanych (H11 i H22 są mniejsze od zera).*)

$Assumptions=Element[{H11,H12,H22,S12},Reals]&&H11≤0&&H22≤0;

H=(({{H11, H12}, {H12, H22}}));

S=(({{1, S12}, {S12, 1}}));

c=(({{c1}, {c2}}));

(*Współczynniki liniowe są określone z dokładnością do stałej, zatem tak naprawdę obchodzi nas ich stosunek (cratio=c1/c2).*)

c=(({{cratio}, {1}}));

(*Oczywiście możemy rozwiązać ten problem na kilka sposobów. Na kartce najłatwiej dla tak małej bazy rozwiązać to z warunku zerowania się wyznacznika.*)

Energie=FullSimplify[Solve[Det[H-En*S]==0,En]]

Wspolczynniki=Table[FullSimplify[Solve[(H-En*S).c==0/.Energie[[i]][[1]],cratio]],{i,1,2}]

(*Można też zastosować metodę, którą wykorzystywaliśmy w zadaniu z metodą Ritza. Ze względu na to, że obliczenia są teraz symboliczne, nie numeryczne, jest to potwornie mało wydajne (Mathematica ciągnie wszystkie symbole zmiennych, pierwiastki, ułamki itp.). Ale "for purely academic purposes..."*)

DiagM[M_]:={{Cos[alfa[M]],Sin[alfa[M]]},{-Sin[alfa[M]],Cos[alfa[M]]}}
alfa[M_]:=1/2ArcCot[(M[[2]][[2]]-M[[1]][[1]])/(2*M[[1]][[2]])]

U=DiagM[S];

SDiag=Simplify[Inverse[U].S.U];

SDiag12=Sqrt[SDiag];

SDiagMinus12=Inverse[SDiag12];

Sm12=Simplify[U.SDiagMinus12.Inverse[U]];

HNew=Simplify[Sm12.H.Sm12];

U=Simplify[DiagM[HNew]];

HDiag=Inverse[U].HNew.U;

(*Te wszystkie ceregiele z linii poniżej są konieczne w celu uproszczenia wzoru (choć nie przeczę, że da się to może zrobić jakoś bardziej elegancko) i uzgodnienia go z postacią otrzymaną powyżej. Diagonalizacja trochę "rozbabrała" końcowy wzór, co nie ma oczywiście miejsca, gdy rozwiązujemy problem numerycznie i wszystkie cząstkowe wyniki są liczbami.*)

Energie=Table[Together[Collect[Together[Apart[TrigReduce[HDiag[[i]][[i]]]]/.Sqrt[a_]->FN[a]],FN[a_]]/.{H11-H22->dH,-H11+H22->-dH}/.dH*FN[a_]:>FN[a*dH^2]/.FN[a_]*FN[b_]:>FN[FullSimplify[Together[a*b]/.dH->H11-H22]]]/.FN[a_]->Sqrt[a],{i,1,2}]

(*Żeby było wydajnie i w zwartym zapisie, można też użyć funkcji wbudowanych w Mathematicę.*)

FullSimplify[Eigensystem[{H,S}]]

(*Przetestujmy teraz szczególne przypadki:*)

(*Degeneracja: H22=H11, ortogonalność stanów: S12=0. *)
ES=Eigensystem[{H,S}/.{H22->H11,S12->0}]

Normalizacja=Table[1/Sqrt[(S[[1]][[1]]^2+S[[2]][[2]]^2+ES[[2]][[i]][[1]]*2S[[1]][[2]])/.{H22->H11,S12->0}],{i,1,2}]

Funkcje=Table[Normalizacja[[i]]*(Psi1+ES[[2]][[i]][[1]]*Psi2),{i,1,2}]

(*Degeneracja, ale brak ortogonalności.*)

ES=Eigensystem[{H,S}/.{H22->H11}]

Normalizacja=Table[1/Sqrt[(S[[1]][[1]]^2+S[[2]][[2]]^2+ES[[2]][[i]][[1]]*2S[[1]][[2]])/.{H22->H11}],{i,1,2}]

Funkcje=Table[Normalizacja[[i]]*(Psi1+ES[[2]][[i]][[1]]*Psi2),{i,1,2}]

(*Brak degeneracji, ale ortogonalność. *)

(*W celu uproszczenia wzorów nazwałem "SQRT" następujący pierwiastek:*)
Sqrt[H11^2+4*H12^2-2*H11*H22+H22^2]

ES=Eigensystem[{H,S}/.{S12->0}]/.Sqrt[H11^2+4*H12^2-2*H11*H22+H22^2]->SQRT

Normalizacja=Table[1/Sqrt[(S[[1]][[1]]^2+S[[2]][[2]]^2+ES[[2]][[i]][[1]]*2S[[1]][[2]])/.{S12->0}],{i,1,2}]

Funkcje=Table[Normalizacja[[i]]*(Psi1+ES[[2]][[i]][[1]]*Psi2),{i,1,2}]