Plot particle motion in potential

1

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

Casimir

Posted 2016-12-18T14:07:20.080

Reputation: 1 251

@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. – Casimir – 2016-12-18T14:20:20.827

Just copy and paste you code from the notebook, select it and hit the {} button. – corey979 – 2016-12-18T14:23:59.293

Exactly what I did. This isn't my first question but I never encountered this problem before. – Casimir – 2016-12-18T14:24:42.193

Can you post the whole question and ignore the message? Or it won't let you post, except for this fragment? – Michael E2 – 2016-12-18T15:13:42.743

@MichaelE2 Correct, it won't let me post more than the above. – Casimir – 2016-12-18T20:30:22.773

I posted a question on meta about this problem. I'll complete this question as soon as the issue can be resolved.

– Casimir – 2016-12-18T20:38:00.490

1It seems you were able to post it. What changed? – Szabolcs – 2016-12-18T20:49:35.993

@Szabolcs I have no idea. As far as I can tell, nothing changed. I copied the markdown text from the meta question back here and it went through without any problems. – Casimir – 2016-12-18T20:51:39.587

Answers

4

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

anderstood

Posted 2016-12-18T14:07:20.080

Reputation: 12 994

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? – Casimir – 2016-12-19T14:10:13.060

@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...). – anderstood – 2016-12-19T15:20:59.587

I don't want to assign values globally, only inside a certain cell (or range of cells). – Casimir – 2016-12-19T15:22:43.257

1@Casimir Then you can use Block[{k=1,m=3}, ...]. – anderstood – 2016-12-19T15:23:17.027

I don't think you can call this particular surface equipotential – BlacKow – 2016-12-19T18:32:24.887

@BlacKow As I understand, $V$ is the total energy which is conserved here (you can see that the trajectory is embedded in the surface). So $V$ corresponds to a isosurface of the total energy. Do you prefer this terminology? – anderstood – 2016-12-19T18:40:38.923

1@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. – BlacKow – 2016-12-19T18:56:03.997

1@BlacKow You're right, edited. – anderstood – 2016-12-19T19:10:27.853

2

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

BlacKow

Posted 2016-12-18T14:07:20.080

Reputation: 6 158