From the Simple Pendulum to Planetary Orbits
We’ve seen in the previous spreadsheet that using the leapfrog technique dramatically improves the accuracy of the spreadsheet in predicting how the pendulum oscillates as time goes on. In this section, we shall extend that technique to the movement of a point mass in two dimensional space. This will enable us to plot the trajectories of planets in the gravitational field of the sun.
The One-Dimensional Simple Harmonic Oscillator
As a simple opening exercise, we consider the case of a one-dimensional simple harmonic oscillator. This is almost the same spreadsheet as the simple pendulum, so open that spreadsheet and save it as 1D_SHO. The simple harmonic oscillator is best thought of as a mass on a spring, the other end of the spring being fixed, and the mass oscillating back and forth on a smooth frictionless table. So starting from the simple pendulum spreadsheet, replace theta by x, replace omega by v, and acceleration by a. Furthermore, replace the acceleration formula, -(g/L)sin theta, by –kx (we’ll take the mass to be one).
This means putting in a fresh set of names: k, x_init, v_init, delta_t, and naming cells in the standard way.
We’ll make one slight change in the formulas: in B19, put x_init, following the previous practice, but in C19 put =v_init+0.5*delta_t*D19 instead of just =v_init. In other words, we’ll put in the halfway jump in the very first cell. In D19, of course, we have just =-k*B19, in other words, -kx, the force (we’re taking unit mass).
On the next row, A20 is just =A19 + delta_t, B20 is =B19+C19*delta_t, C20 is =C19+D20*delta_t and D20 is =-k*B20.
We are now ready to roll, and copy A20 through D20 down 200 rows, then create a graph showing the position of the oscillator as a function of time. This is less interesting than the simple pendulum. It is always a sine curve, independent of amplitude. Of course, that isn’t totally physically realistic, any spring will ultimately break.
Save this spreadsheet.
The Two-Dimensional Harmonic Oscillator
Physically, a two-dimensional harmonic oscillator is just a point mass constrained to move in a plane, with a force directed towards the origin, and of magnitude directly proportional to r, the distance from the origin.
Representing the force kr directed towards the origin by a vector, we see that this can be written as the sum of two vectors, one parallel to the x axis and one parallel to the y axis, of lengths kx and ky respectively, and both directed towards the origin. We can therefore write the vector equation of motion of the point mass, F = ma, in terms of its components, Fx = max, Fy = may. For unit mass this gives ax = -kx, ay = -ky.
This means that the x and y motions are independent one-dimensional oscillators! At this point, we open up the 1D_SHO spreadsheet and re-save it as 2D_SHO. We’re just going to add in extra columns for the motion in the y direction.
First we need extra variables y_init and v_y_init, and we should rename v_init as v_x_init. So put these new names in the system.
Now insert new columns for the variables y and v_y, and don’t forget to relabel the v column as v_x.
By copying the formulas over from the x columns to the corresponding y columns, and making appropriate adjustments (replace v_x_init by v_y_init) a 2D spreadsheet can be constructed.
Plot x against y in a scatter plot. In general you will see an ellipse, centered at the origin. If you take delta_t very large, it will break up, signalling the limitations of the numerical method.
Actually, the two-dimensional oscillator described above is in a sense rather trivial, in that it is really just the sum of two one-dimensional oscillators, the movements in the x and y directions are completely independent. However, it is well worth constructing the spreadsheet, because it is easily adapted to much more interesting systems, such as the movements of planets under the inverse-square gravitational force of the sun. This force we will write for simplicity as k/r2 directed towards the center, so k = GMm.
Since we will work with accelerations in the x and y directions, we must write this vector force in terms of its x and y direction components. Draw a diagram to convince yourself that these components have magnitudes (x/r)k/r2 and (y/r)k/r2, or kx/r3 and ky/r3, where of course r2 = x2 + y2.
All we need do to make out 2D_SHO spreadsheet describe instead planetary orbits is to replace the formulas –kx and –ky for accelerations with –kx/r3 and –ky/r3.
So, open the 2D_SHO spreadsheet and re-save it as Planet. Now, in F19 write:
and in G19 write
Copy these formulas down 200 rows.
Adjust the values of the initial parameters, including k and delta_t, until you get an elliptical planetary orbit. Unlike the two-dimensional simple harmonic oscillator, this will not be centered at the origin.