Copyright © Michael Richmond.
This work is licensed under a Creative Commons License.
Our problem -- a cannonball shot towards the Moon -- illustrates one of the big challenges for numerical methods: dealing with forces and accelerations which vary over wide ranges.
When the projectile is flying through the Earth's atmosphere at high speed, it encounters very large forces of air resistance. For example, if we use this simple approximation
2
F(air) = C * (density) * (area) * (speed)
and some reasonable parameters, we end up with
2
F(air) = 0.02 * (1.21 kg/m^3) * (12 m^2) * (10,000 m/s)
= 29,000,000 Newtons
Compare that to the force of gravity on our projectile, which has a mass of 10 tonnes:
F(grav) = mass * g
= 98,000 Newtons
Yow! The force of air resistance overwhelms the gravitational force.
On the other hand, consider the force on the same projectile once it has gone above the Earth's atmosphere. At an altitude of 200 km, for example, air resistance is negligible, and the force of gravity drops to
F(grav) = G * M(earth) * m(ship) / ( total R )^2
= 92,400 Newtons
When the cannonball is halfway to the Moon, this force is very much weaker still:
F(grav) = G * M(earth) * m(ship) / ( 1.92E8 m )^2
= 108 Newtons
The total force on the projectile varies by many orders of magnitude (by a factor of roughly 270,000, to be specific). What does this mean?
If we use a fixed timestep for the entire simulation, we will have to pick a small one, like dt = 0.01 seconds, to avoid errors during the first few seconds of flight. But that's a pain in the neck:
Q: If we use a stepsize of 0.01 seconds,
and the total flight takes 5 days,
how many steps must we make?
One way to avoid such long execution times is to vary the size of the timestep: use a small step when the accelerations are large, and a large step when the accelerations are small. There are two questions we must answer:
The basic idea is to use the size of the acceleration of the cannonball as an indicator: big accelerations mean we need small timesteps. Actually, that's not quite true: if we use Heun's method, then we can handle even a large acceleration AS LONG AS IT IS CONSTANT. The trouble occurs if the acceleration changes quickly and significantly from one timestep to the next.
Fortunately, Heun's method provides a very simple way to check signficant changes: it includes the prediction of the acceleration we'll experience in a future timestep. So, we can simply compute a fractional change:
( current accel - predicted accel )
frac_change = abs ( ------------------------------- )
( current accel
If the fractional change in acceleration is larger than some threshold, we might have to choose small timesteps; if it is smaller than some threshold, we might try using larger timesteps.
A simple algorithm for varying the timestep might look sort of like this:
compute frac_change in acceleration
if (frac_change > thresh_1)
new_timestep = 0.5 * timestep;
if (frac_change < thresh_2)
new_timestep = 2.0 * timestep;
otherwise, leave timestep as-is
Note that this algorithm allows the timestep to float into any range; in some circumstances, you might want to set minimum and maximum permitted values.
What are the appropriate threshold values? Good question. They will probably depend exactly on the nature of the calculations you are performing. The next section may help you to pick reasonable values.
If you know the exact solution to your problem, then you can compare the simulated results to the exact solution and decide if the errors are acceptable. In this week's case, however -- a projectile flying through the atmosphere -- we don't know the exact solution. Rats.
Could we use conserved quantities, such as total energy or angular momentum, as a check? In some situations, yes -- but in this situation, no: the total energy drops a bit as the projectile suffers air resistance.
Here are two things we can do:
Of course, just about any technique will yield some small change in energy. You must decide, for each particular situation, what sort of change is important, and what sort of change is not important. For example, you might decide that an error which changes the travel time by a few minutes -- out of a total of several days -- is not important. Once you've made this decision, you can start fiddling with the thresholds and limits on your timestep to make sure that your final result falls within the acceptable range.
Here's a version of the program for last week's assignment: it simulates the motion of a cannonball without air resistance, but including a changing gravitational force.
You can call it like so:
cannonball_a(initial_speed, timestep, "heun", "test.out")
This will write diagnostic messages to the screen as the program runs, and also write a more compact set of output to the file "test.out".
Here's a slightly different program: it includes air resistance, but NOT the changing gravitational force from the Earth.
You can call it like so:
cannonball_b(initial_speed, timestep, "heun", "test.out")
Copyright © Michael Richmond.
This work is licensed under a Creative Commons License.