# Using Loop to find a specific Energy

I am given a Hamiltonian:

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


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

H = 9.28*10^-6


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

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


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

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


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

• I am a little confused by your notation. Are you perhaps looking for: NSolve[H == 9.28*10^-6, {c, c, c', c'}, Reals] ? It gives c -> 0.727438, c -> -4.12957, c' -> 3.98141, c' -> 1.89354. Apr 26, 2022 at 15:42
• @Domen NSolve gives only one solution. As there are infinitely many solutions for this problem, I was looking for a finite set of solutions(in the desired range) by using a loop. Apr 27, 2022 at 5:46

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

vars = {c', c', c, c};
FindInstance[H == 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 == 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.

• Thank you for enlightening me with the "FindInstance" command :) I needed these initial conditions for a Poincare section problem, that I have been stuck with for weeks: mathematica.stackexchange.com/questions/265876/… I will try to use these and see if I get any success. Thanks again! Apr 27, 2022 at 15:31

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

(* {{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}} *)

• Thank you for taking out some time from your busy schedule to solve this problem. It is really helpful. Actually, I wanted to find these random solutions for a Poincare section problem, I have been stuck in for weeks: mathematica.stackexchange.com/questions/265876/… I wished to find the poincare section for the given fixed energy in order to reproduce results from this paper(page 5): arxiv.org/pdf/1803.06756.pdf . I will try to use your algorithm and see if I get any success. Thank you again :) Apr 27, 2022 at 14:07