Need angular momentum error < 1% - Don't know why error values aren't changing w/ changing time step

Hello all,

Before I jump in to the details here's my code -- https://pastebin.com/y6wcUHg6

I'm working on some code that simulates an orbit and then calculates the kinetic, potential and total energy of the orbit and displays them. This was originally skeleton code for a class in which we needed to update from the Euler Method to the 2nd order Runge Kutta method and also calculate the angular momentum as well as the energies previously mentioned and plot that as well.

----------------------------------------------------------------------------------------------------------------------------------------------------------------

We know in an ideal earth orbit (without the influence of the moon) the angular momentum and total energy are conserved quantities, thus the initial values should equal the final values. We are trying to find a time step size for which the fractional change in these quantities (discussed below) is less than 1%. Now I am tasked with finding how small a step size I is needed in order to find a less than 1% change in the energy for a typical earth orbit -- using r0 = 1 AU. v0 = 2*pi AU/Yr and then a combination of nStep (# of steps) and tau (time step size) such that the product of the 2 always equals 1 (i.e., nStep = 100 and tau = .01). And then finding how small a tau I need in order to find a less than 1% change in the angular momentum as well.

----------------------------------------------------------------------------------------------------------------------------------------------------------------

In my code I haven't yet written anything to display the fractional change in these quantites {something like an fprintf statement displaying, say, 100*abs(L(final)-L(initial)/L(initial))} but upon inspection I start with nStep = 100 and tau = .01 and get a potential energy around -40 a kinetic around +20 and total around -20. These energies are technically not showing any error on them as they are constant in time but we did a similar exercise recently and that was not what we observed, we observed the total energy change growing closer and closer to zero as the time step decreased. I can provide that code as well, if necessary at all. Then I get a 1x100 matrix of angular momentum values which has an initial value of around 31 and a final value around 3.2*10^5. The graph I get reflects this.

When I decrease the time step by factors of 10 (and consequently increase nStep by the same factors of 10) I continue to get the same values for the energies and for every factor of ten I decrease the time step by the final value for the angular momentum seemingly increases by a factor of 1000.

I don't believe that increasing a time step would give me a fractional angular momentum change closer to 0 and I'm not sure why my energies aren't varying in their initial and final values (or any for that matter) and their fractional energy change converging to zero as the time steps approach infinity. Any help is greatly appreciated!! Have a lovely evening :)

[Kyle O'Connell]

I've been doing some fiddling around and I'm seeing that the only thing that changes the change in the initial and final values is changing how long the for loop runs for (i.e., changing nStep).. but this seems trivial as it's basically saying you have this iteration which takes an initial value then applies some formula to it for a number of steps, each step ostensibly moving the current value farther away from the initial, thus running for less time afford for less "drift" from the initial value. Could it be that something this simple is being asked and I simply need to find a combination of nStep and Tau which allows for a short enough time such that the change in these quanitites is less than 1%? Or am I missing something fundamental about the question (or orbits)? Thanks!

Answer 1

[David Goodmanson]

Hello Kyle,

Right now you are calculating the orbit using r and v and you don't need any other variables to do that. However, if you wanted to do a parallel calculation for L, you would have at each time step

L = L + tau*torquemid

similar to what you are doing with r and v.

If there were a torque exerted on the planet as it moves around the sun, then the angular momentum would be changing with time. However, since the force is a central force in the direction of -r, the torque (r cross F) is zero. You expect L to be a constant. So the array 'angular' should just be a check on how well the orbit calculation is going. There are three lines of code involving w0,L and angular, and if you replace that stuff with

angular(iStep) = mass*det([rmid; vmid]);

then you can see how well your code is doing in preserving L.