I am trying to bootstrap the harmonic oscillator in Mathematica using its built in SDP solver. I have the following code
k = 55;
reqRel[s_] :=
2*x*ToExpression["x" <> ToString[s - 2]]*(s - 1) +
1/4*(s - 1) (s - 2) (s - 3) ToExpression["x" <> ToString[s - 4]] -
s ToExpression["x" <> ToString[s]] == 0;
rels = Table[reqRel[2 j], {j, 1, k}];
For[i = 0, i <= Length[rels], i++,
rules =
If[i == 0, { ToExpression["x" <> ToString[0]] -> 1},
FullSimplify[
Join[rules,
First[Solve[rels[[i]] /. rules,
ToExpression["x" <> ToString[2 i]]]]]]]
];
(*This section implements the moment recursion relation for the hamiltonian as a rule to be used later*)
matt = Table[
If[OddQ[i + j], 0,
If[i == 0 && j == 0, 1,
ToExpression["x" <> ToString[i + j]]]], {i, 0, k}, {j, 0, k}
];
(*Setting up the matrix with even moments of x as its even numbered (even sum of row and columns)*)
funcmatt[e_] := (matt /. rules) //. {x -> e}
(*setting up the matrix as a function to be used inside the SDP*)
iterationAtk[e_] := SemidefiniteOptimization[
-t,
VectorGreaterEqual[{funcmatt[e] - t IdentityMatrix[k + 1],
0}, {"SemidefiniteCone", k + 1}],
{t \[Element] Reals}, "PrimalMinimumValue"]
data = {};
acceptedPoints = {};
Do[
points = iterationAtk[e];
AppendTo[data, {e, points}];
If[points >= 0, AppendTo[acceptedPoints, {e, points}]];
, {e, 0, 10, 0.1}]
acceptedPoints
data
(*I iterate through energies, since the constraint doesn't have the proper curvature with respect to the energy*)
I feel like I am missing something here conceptually, most likely in how I am implementing the SDP. I also run a test by just checking the eigenvalues of the matrix and doing the same itteration:
datanew = {};
Do[
If[Min[Eigenvalues[funcmatt[e]]] >= 0, AppendTo[datanew, e];
Print[{Min[Eigenvalues[funcmatt[e]]], e}]],
{e, 0, 10, 0.1}
]
datanew
In both cases clearly allowed energies like 0.5 and 1.5 are not allowed according to the code. (I am working in natural units, so hbar=m=omega=1). It is worth mentioning that I am in general following this publication.If needed, I have also uploaded the notebook here(github).Any help is appreciated.