0
$\begingroup$

Consider the anisotropic harmonic potential in two dimensions $(q_1,q_2)$ given by

$$ V(q_1,q_2) = \frac{m}{2} \, q_1^2 + \frac{k}{2} \, q_2^2, $$

or V = m/2 q1^2 + k/2 q2^2; in Mathematica.

The Newtonian e.o.m.s of a particle moving through this potential are

$$ \ddot{q}_1 = -q_1, \qquad \ddot{q}_2 = -\omega^2 \, q_2, $$

where $\omega = \sqrt{k/m}$ is the angular frequency of oscillations in the $q_2$-direction. Given the initial conditions $q_i(0) = q_{i,0}$ and $p_i(0) = p_{i,0}$, the e.o.m.s are solved by

$$ \begin{aligned} q_1(t) &= q_{1,i} \cos(t) + \frac{p_{1,i}}{m} \, \sin(t),\\ q_2(t) &= q_{2,i} \cos(\omega t) + \frac{p_{2,i}}{m \omega} \, \sin(\omega t), \end{aligned} \qquad \text{with $p_i = m \dot{q}_i$.} $$

In Mathematica:

DSolve[{q1''[t] == -q1[t], q2''[t] == -ω^2 q2[t], q1[0] == q10,
q1'[0] == p10/m, q2[0] == q20, q2'[0] == p20/m}, {q1[t], q2[t]}, t]
//FullSimplify

{{q1[t] -> q10 Cos[t] + (p10 Sin[t])/m, q2[t] -> q20 Cos[t ω] + (p20 Sin[t ω])/(m ω)}}

A 3d plot of the potential looks like this.

Plot3D[V /. {m -> 1, k -> 3}, {q1, -5, 5}, {q2, -5, 5},RegionFunction -> Function[{q1, q2}, m/2 q1^2 + k/2 q2^2 <= 12 /. {m -> 1, k -> 3}]]

enter image description here

What I would like to do now is draw the particle as a ball moving through this potential with a fixed energy. Any help would be much appreciated.


Some remarks about the physics behind this simulation: The particle always reaches a certain height both in $q_1$- and $q_2$-direction before rolling back down again and up the other side. This is because the Hamiltonian $H = \frac{p_1^2}{2 m} + \frac{p_2^2}{2 m} + V(q_1,q_2)$ does not couple the degrees of freedom in $q_1$- and $q_2$-direction. Therefore, the total energies $E_1$ and $E_2$ available in dimensions $q_1$ and $q_2$ are conserved separately.

For long times $t \to \infty$ and irrational angular frequency $\omega \notin \mathbb{Q}$ (which ensures that the trajectory never closes, thus making the system ergodic), the particle's trajectory should therefore trace out a rectangle $R$ whose length and width are determined by $E_1$ and $E_2$.


Update: With anderstood's and BlacKow's help, I was able to piece together this solution that does exactly what I want.

V = m/2 q1^2 + k/2 q2^2; \[Omega] = Sqrt[k/m];
sol = {q1[t], q2[t], m/2 q1[t]^2 + k/2 q2[t]^2} /. 
   DSolve[{q1''[t] == -q1[t], q2''[t] == -\[Omega]^2 q2[t], 
     q1[0] == q10, q1'[0] == p10/m, q2[0] == q20, 
     q2'[0] == p20/m}, {q1[t], q2[t]}, t] // FullSimplify

Manipulate[Block[{m = 1, k = 5, q10 = 5, p10 = 1, q20 = 2, p20 = 1},
  surf = Plot3D[V, {q1, -7, 7}, {q2, -5, 5}, 
    RegionFunction -> Function[{q1, q2}, m/2 q1^2 + k/2 q2^2 <= 25], PlotStyle -> Opacity[0.5]];
  Show[surf, Graphics3D@{Blue, Ball[sol /. {t -> tf}, 0.4]}, 
   ParametricPlot3D[sol, {t, 0, tf}]]], {tf, 0.1, 100}]

enter image description here

$\endgroup$
  • $\begingroup$ @corey979 Sorry, I'm having trouble posting this question. I keep getting the error Your post appears to contain code that is not properly formatted as code. Please indent all code by 4 spaces using the code toolbar button or the CTRL+K keyboard shortcut. For more editing help, click the [?] toolbar icon. Looking on Meta, it appears this issue has come up before. $\endgroup$ – Casimir Dec 18 '16 at 14:20
  • $\begingroup$ Just copy and paste you code from the notebook, select it and hit the {} button. $\endgroup$ – corey979 Dec 18 '16 at 14:23
  • $\begingroup$ Exactly what I did. This isn't my first question but I never encountered this problem before. $\endgroup$ – Casimir Dec 18 '16 at 14:24
  • $\begingroup$ Can you post the whole question and ignore the message? Or it won't let you post, except for this fragment? $\endgroup$ – Michael E2 Dec 18 '16 at 15:13
  • 1
    $\begingroup$ It seems you were able to post it. What changed? $\endgroup$ – Szabolcs Dec 18 '16 at 20:49
4
$\begingroup$

First, store the 3D trajectory (in the space $(q_1,q_2,V)$):

sol = {q1[t], q2[t], m/2 q1[t]^2 + k/2 q2[t]^2} /. 
   DSolve[{q1''[t] == -q1[t], q2''[t] == -\[Omega]^2 q2[t], 
     q1[0] == q10, q1'[0] == p10/m, q2[0] == q20, 
     q2'[0] == p20/m}, {q1[t], q2[t]}, t] // FullSimplify

Plot the potential surface:

surf = Plot3D[
  V /. {m -> 1, k -> 3, \[Omega] -> Sqrt[k/m]}, {q1, -5, 5}, {q2, -5, 
   5}, RegionFunction -> 
   Function[{q1, q2}, m/2 q1^2 + k/2 q2^2 <= 12 /. {m -> 1, k -> 3}]]

Then select some initial conditions and plot the trajectory using ParametricPlot3D. Show it together with the surface using... Show:

p10 = 3; p20 = 1; q20 = -1.5; q10 = 2;
traj = ParametricPlot3D[
   sol /. \[Omega] -> Sqrt[k/m] /. {m -> 1, k -> 3}, {t, 0, 50}, 
   PlotStyle -> Red];
Show[pot, traj]

enter image description here

$\endgroup$
  • $\begingroup$ Looks good! One question: Is there a way to avoid repeating the replacement /. {m -> 1, k -> 3} everywhere m or k come up? For example, could they be assigned specific values all throughout a cell? $\endgroup$ – Casimir Dec 19 '16 at 14:10
  • $\begingroup$ @Casimir I you want to keep these values throughout the notebook, just add at the beginning: m=1; k=3;. And you should also avoid using "redundant" variables, i.e. $\omega$, which depends on k and m (I guess...). $\endgroup$ – anderstood Dec 19 '16 at 15:20
  • 1
    $\begingroup$ @Casimir Then you can use Block[{k=1,m=3}, ...]. $\endgroup$ – anderstood Dec 19 '16 at 15:23
  • 1
    $\begingroup$ @anderstood equipotential surface is a surface where potential energy is the same, so it's defined by $V(q) = const$; it will be an ellipse in $(q_1,q_2)$ coordinates. And here $V$ is not a total energy but, potential energy only. $\endgroup$ – BlacKow Dec 19 '16 at 18:56
  • 1
    $\begingroup$ @BlacKow You're right, edited. $\endgroup$ – anderstood Dec 19 '16 at 19:10
2
$\begingroup$

Using @anderstood answer we can play with graphics to make the surface transparent and plot ball movement:

V[q1_, q2_] := m/2 q1 q1 + k/2 q2 q2
surf = Plot3D[
  V[q1, q2] /. {m -> 1, k -> 3, \[Omega] -> Sqrt[k/m]}, {q1, -5, 
   5}, {q2, -5, 5}, 
  RegionFunction -> 
   Function[{q1, q2}, m/2 q1^2 + k/2 q2^2 <= 12 /. {m -> 1, k -> 3}], 
  Mesh -> None, 
  ColorFunction -> 
   Function[{z}, Opacity[0.4, #] &@ColorData["TemperatureMap"][z]]]

traj[t_] := 
 Evaluate[Flatten@sol /. \[Omega] -> Sqrt[k/m] /. {m -> 1, 
     k -> 3} /. {p10 -> 3, p20 -> 1, q20 -> -1.5, q10 -> 2}]
frames = Table[
   Show[surf, Graphics3D@{Red, Ball[traj[t], 0.2]}], {t, 0, 10, 
    0.1}];
Export["Documents/animBall.gif", frames]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.