Elliptical Discussions
It’s the ninth week of Semester 2 and I’m coming to the end of lectures and laboratory sessions for my Computational Physics module; for the remaining three weeks (plus the Easter vacation) the students in this class will be mainly concentrating on the mini-projects that form part of the assessment.
This afternoon, though we had a session on how to transform higher-order differential equations into sets of coupled first-order ODEs suitable for vectorization and consequent solution using standard techniques. The problem we focussed on today was the simple problem of orbital motion of a test particle under the gravitational force in plane-polar coordinates, which can be prepared for physical solution thusly:
This sort of thing reminds me of my undergraduate theory project at Cambridge, where I did a similar thing to solve the equations governing the action of a four-level laser, though that was in Fortran rather than Python. In my own solution I used Python’s off-the-shelf solver odeint.
I like the orbital motion problem a lot because it’s a bit more than a coding exercise: students have to think about how to choose initial data, how to test the their code and interpret the results. Even before that there’s the issue of what units to use; SI units are a bit daft for astronomical problems. For solar system calculations it makes sense to use Astronomical Units for distances and years for time; in such a system it’s easy to work out that GM is just 4π2, which avoids having to deal with ridiculously large or ridiculously small numbers.
Anyway, the fun thing about this lab was that once everyone had got their code working they could try setting initial data to get a circular orbit as a special case, explore how the shape of elliptical orbits depends on the input data, how to make an unbound orbit, and so on. It’s important to understand the output of a numerical calculation in terms of basic physical principles. All that led to a discussion in class of solar system exploration, transfer orbits, what would happen if the mass of the Sun suddenly changed, or if G was a function of time, and lots of other things.
I find sessions like this that encourage students to explore problems themselves very rewarding and I think they add a valuable extra dimension to standard teaching formats. I hope the projects that they’ll be doing from now on – involving topics in areas ranging from atomic physics, cosmology, particle physics and climate science, and done in groups – will provoke even more discussion of this type.
Follow @telescoper
In the Dark 
April 8, 2022 at 10:22 am
I think your approach to units here leaves something to be desired. When I have students doing projects with any kind of numerical equation solving, I emphasise to them that computers have no idea about units and can only deal with pure numbers. So the first step is always to make the equations explicitly dimensionless before you set about solving them. In this case, I would probably define R = r / r_0, where r_0 is the starting value. Then r_0^3/GM = T^2 is the square of a time unit. So if we divide through by this and work with t / T as a time coordinate, then the equations to be solved have in effect GM = 1. But I think it’s better conceptually to show explicitly the construction of dimensionless variables rather than to talk about “the units of GM”.
April 8, 2022 at 11:16 am
Reasonable point. In fact I did this with another problem – the comparison between the SHM solution for small angles and the full pendulum solution with sin theta instead of theta by using sqrt(l/g) to make the time coordinate dimensionsless.
In the orbits case I chose to do it the way I did primarily because it was very straightforward to choose initial data that would give a circular orbit at the Earth’s radius and check the period came out right as T=1 for a=1. The initial velocity is just 2pi in this system.