Haven’t updated my blog in a while.

Basically having a bit more fun with VPython. I’ve been simulating solar system, pendulum and spring-mass system.

My usual approach is to to consider all the forces, express them as vector and sum them up. Then calculate the acceleration, then the velocity, and from that the position.

For example, in the two-body solar system where the sun and planet pull on each other, it is relatively straightforward to calculate the gravitational pull each body exerts on the other, then calculate acceleration, velocity and update the positions.

The VPython interpreter I had downloaded, however, comes with a number of examples, including a two-body problem which I couldn’t understand at first, since it was compact and made use of momentum rather than forces. I was sufficiently rusty in first year physics that I couldn’t figure it out.

When I wrote my own two-body solar system, I used my basic approach of calculating forces, acceleration and velocity and position. I found that the system tended to wander off in one direction. Then I remembered from the example program the requirement that the sun’s momentum be equal and opposite to the planet’s momentum. That is, the sum of momentum should be zero. I thought about that for a while. I realized that if the sum of momentum was nonzero, the whole system would have excess momentum and move in a direction. I had been wondering about that requirement.

When I carefully input the initial velocities of the sun and planets so that their momenta would sum to zero, the solar system stopped drifting.

Then I thought about the two-body system some more. I had been teaching undergrads for a while and it refreshed me on the basics. I recalled that force is the rate of change of momentum.

*dF = dp/dt*

So *dp = F dt*

p – *p*_{0} = -G Mm **r **/ r3 dt

thus p = *p*_{0} -G Mm **r **/ r3 dt

Momentum of sun = -Momentum of planet, so I only had to do the calculation once.

As for position, *ds _{planet} = p/m_{planet} dt *and

*ds*

_{sun}= p/m_{sun}dtSince ds = s_{new} – s_{old} , we get *s _{new} = s_{old} + p/m dt .*

The result was a very compact program that did the necessary calculation in one line, and updated the positions of the bodies in two more lines.

Damn it was elegant. I’m sorry to say my proficiency still isn’t at the level of coming up with such a compact and elegant algorithm, I tend to take more steps than I ought to… although I still get there at the end.