Using Loop to find a specific Energy

Question

I am given a Hamiltonian:

H[t_] = c[0]'[t]^2 + c[1]'[t]^2 + 4.66 c[0][t] c[0]'[t]^2 + 
   1.32 c[0][t] c[1]'[t]^2 - 7.57 c[0]'[t] c[1][t] c[1]'[t] - 
   1.40 c[0][t]^2 + 7.57 c[1][t]^2 - 7.11 c[0][t]^3 - 
   35.3 c[0][t] c[1][t]^2;

I wish to find a suitable set of values for the different variables such that I get the answer as:

H[0] = 9.28*10^-6

In order to go around solving this problem, I tried using the RandomReal function of Mathematica:

c[0][0] = RandomReal[{-0.01, 0.01}]; 
Derivative[1][c[0]][0] = RandomReal[{-0.01, 0.01}]; 
c[1][0] = RandomReal[{-0.01, 0.01}]; 
Derivative[1][c[1]][0] = RandomReal[{-0.01, 0.01}];

And then simply use the While loop to get the desired results:

While[H[0] = 9.28*10^-6, Print[c[0][0]]; 
 Print[Derivative[1][c[0]][0]]; Print[c[1][0]]; 
 Print[Derivative[1][c[1]][0]]]

but the code doesn't work. How shall I rectify my error?

Answer 1

If you just want a set of solutions that satisfy your initial condition, you could use FindInstance:

vars = {c[0]'[0], c[1]'[0], c[0][0], c[1][0]}; 
FindInstance[H[0] == 9.28*10^-6, vars, Reals, 5]

The number 5 can be adjusted to yield different numbers of solutions. If you want the results to be constrained to a particular range, you can do that as well; it just takes longer:

FindInstance[{H[0] == 9.28*10^-6 && And @@ Thread[-0.01 < vars < 0.01]}, vars, Reals, 5]

Note that these results are not "random", in the sense that Mathematica will return the same results each time you run the code. They are also not distributed according to any particular probability distribution, but are rather determined by whatever black box the good people at Wolfram use inside of this code.

Answer 2

Let's first fix your approach. I have slightly changed the naming of the variables, but the main difference is in using SetDelayed (:=) instead of Set (=). With this, everytime you evaluate the variable, it will get a new, random value.

Clear[h, h0, c0, c1, dc0, dc1];

h = dc0^2 + dc1^2 + 4.66 c0 dc0^2 + 1.32 c0 dc1^2 - 7.57 dc0 c1 dc1 - 
   1.40 c0^2 + 7.57 c1^2 - 7.11 c0^3 - 35.3 c0 c1^2;
h0 = 9.28*10^-6;
c0 := RandomReal[{-0.01, 0.01}];
dc0 := RandomReal[{-0.01, 0.01}];
c1 := RandomReal[{-0.01, 0.01}];
dc1 := RandomReal[{-0.01, 0.01}];

While[h != h0]

Run this code and wait. You will see that it will not terminate in any reasonable time. This is because it is extremely unlikely that if you randomly choose four numbers, they would obey your non-trivial equation. We need another approach.

Let's first try to solve the equation.

Clear[c0, c1, dc0, dc1];

Solve[Rationalize[h == h0], {c0, c1, dc0, dc1}]

We see that there are several branches of solution with different dimensionalities, but the largest one is 3-dimensional (3 degrees of freedom). Since you haven't given any constraints or requirements about the "randomness" of the solution, we can use this branch: generate randomly three of the variables, calculate the fourth one, and repeat this until the fourth one is in the desired interval. Below is an example of a (rather clumsy) implementation.

Clear[generateRandomSolutions];

generateRandomSolutions[n_ : 1] := Module[{
    h = dc0^2 + dc1^2 + 4.66 c0 dc0^2 + 1.32 c0 dc1^2 - 7.57 dc0 c1 dc1 - 
      1.40 c0^2 + 7.57 c1^2 - 7.11 c0^3 - 35.3 c0 c1^2,
    h0 = 9.28*10^-6,
    interval = {-0.01, 0.01},
    ret = Nothing,
    eq, insideInterval, res, sol,
    rets = {}
   },

  eq = Rationalize[h0 == h];
  insideInterval = IntervalMemberQ[Interval[interval], #] && # \[Element] Reals &;

  Do[
   While[ret == Nothing || !insideInterval[Values@Last@ret],
    sol = First@Solve[eq, dc1];
    res = MapThread[Rule, {{c0, c1, dc0}, RandomReal[interval, 3]}];
    ret = res~Join~(sol /. res);
    ];
   AppendTo[rets, ret];
   ret = Nothing, n];
  
  rets
]

generateRandomSolutions[3]

(* {{c0 -> 0.00614171, c1 -> 0.000418391, dc0 -> -0.00662086, dc1 -> -0.00415999}, 
 {c0 -> 0.00828939, c1 -> -0.00166806, dc0 -> -0.0043044, dc1 -> -0.00829657}, 
 {c0 -> -0.00544832, c1 -> -0.000750933, dc0 -> 0.00313956, dc1 -> -0.00600589}} *)