0
$\begingroup$

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.

$\endgroup$

0

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.