C++
Optional prize problem
TASK DESCRIPTION
If you finish the rest of the lab early and would like to enter for
the Part 1B computing prize then you may attempt this exercise. There
are no marks for this exercise. In order to be considered for the
Prize, your solution must be received by Prof Roberto Cipolla by one week
after the end of the Lent Full Term 2010, in the form specified at the
end of this document.
Numerically integrate the orbits of celestial bodies in order to predict when solar eclipses will occur. Given below are the positions, velocities, masses, and radii of the Moon, Earth, and Sun at midnight on the beginning of 7 January, 2008. The Cartesian coordinates are given relative to the barycentre of the solar system, with axial directions determined by the orbital plane of the earth. Positions are given in astronomical units (AU), and velocities are given in AU per day:
| Body | x (AU) | y (AU) | z (AU) | vx (AU/day) | vy (AU/day) | vz (AU/day) | mass (kg) | radius (m) |
| Moon | ||||||||
| Earth | ||||||||
| Sun |
The formula for gravitational attractive force between bodies A and B is given by
FAB = G⋅mA⋅mB dAB2
where G is the gravitational constant, mX is the mass of body X in kilograms, and dAB is the distance between the two bodies in meters. Use these values for G and the astronomical unit:
| G | = | 6.67300×10-11 m3⋅kg-1⋅s-2 |
| 1 AU | = | 1.49598×1011 m |
RUNGE-KUTTA INTEGRATION
The fourth-order Runge-Kutta method (RK4) for numerical integration evaluates the derivative at 4 different points at each step. Consider the initial value problem:
| y' | = | f(y) |
| y(t0) | = | y0 |
Choose a step size h. Then the RK4 update step is given by:
| k1 | = | f(yn) | ||||
| k2 | = | f(yn + ½⋅h⋅k1) | ||||
| k3 | = | f(yn + ½⋅h⋅k2) | ||||
| k4 | = | f(yn + h⋅k3) | ||||
| yn+1 | = |
|
The per-step error of the RK4 method is on the order of h5. In contrast, the Euler method has a per-step error of h2, requiring much smaller h for the same error bound.
EVALUATION
Try the Euler method, the RK4 method, and perhaps other numerical integration methods (e.g. Gaussian Quadrature) to compute the trajectories of the bodies, and calculate the expected times for solar eclipses. Try varying the step size of integration and the angular alignment tolerance for what should be considered an eclipse. A list of predicted eclipse dates is given by Wikipedia. How far into the future can your program correctly predict eclipses relative to the list? How does this vary with step size and integration method?
ANSWER SUBMISSION
Your answer must be submitted entirely on paper, and sent to Prof Roberto Cipolla through the internal mail (hand it in at the Baker Building Enquiry office). It should contain the following parts:
- Description of the technique you have used to solve the problem.
- Source code for the programs you have written. Please include appropriate comments.
- The estimated eclipse times and any other computed data.
- Comments and conclusions.
Submissions must be received in this form by one week after the end of Lent Full Term 2010 if they are to be considered for the prize.
Back to main 1B C++ page