ODEs

Takin from Newmon’s Chapter 8

8.7 Trajectory with air resistance

Many elementary mechanics problems deal with the physics of objects moving or flying through the air, but they almost always ignore friction and air resistance to make the equations solvable. If we’re using a computer, however, we don’t need solvable equations.

Consider, for instance, a spherical cannonball shot from a cannon standing on level ground. The air resistance on a moving sphere is a force in the opposite direction to the motion with magnitude

F= \frac{1}{2} \pi R^2 \rho C v^2

where R is the sphere’s radius, \rho=1.22 kg/m3 is the density of air, v is the velocity, and C is the so-called coefficient of drag (a property of the shape of the moving object, in this case a sphere). Starting from Newton’s second law, F=ma, show that the equations of motion for the position (x,y) of the cannonball are

\ddot{x} = - {\pi R^2\rho C\over2m}\, \dot{x} \sqrt{\dot{x}^2+\dot{y}^2}
\ddot{y} = - g - {\pi R^2\rho C\over2m}\, \dot{y}\sqrt{\dot{x}^2+\dot{y}^2},

where m is the mass of the cannonball, g is the acceleration due to gravity, and \dot{x} and \ddot{x} are the first and second derivatives of x with respect to time.
Change these two second-order equations into four first-order equations using the methods you have learned, then write a program that solves the equations for a cannonball of mass 1kg and radius 8cm, shot at 30° to the horizontal with initial velocity100m/s. The density of air is \rho=1.22\,\textrm{kg}\,\textrm{m}^{-3} and the coefficient of drag for a sphere is C=0.47. Make a plot of the trajectory of the cannonball (i.e., a graph of y as a function of x).
When one ignores air resistance, the distance traveled by a projectile does not depend on the mass of the projectile. In real life, however, mass certainly does make a difference. Use your program to estimate the total distance traveled (over horizontal ground) by the cannonball above, and then experiment with the program to determine whether the cannonball travels further if it is heavier or lighter. You could, for instance, plot a series of trajectories for cannonballs of different masses, or you could make a graph of distance traveled as a
function of mass. Describe briefly what you discover.

Exercise 8.8 Space garbage

A heavy steel rod and a spherical ball-bearing, discarded by a passing spaceship, are floating in zero gravity and the ball bearing is orbiting around the rod under the effect of its gravitational pull:

rod

For simplicity we’ll assume that the rod is of negligible cross-section and heavy enough that it doesn’t move significantly, and that the ball bearing is orbiting around the rod’s mid-point in a plane perpendicular to the rod.
Treating the rod as a line of mass M and length L and the ball bearing as a point mass m, show that the attractive force F felt by the ball bearing in the direction toward the center of the rod is given by

F = {GMm\over{L}} \sqrt{x^2+y^2} \int_{-L/2}^{L/2} {dz\over{(x^2+y^2+z^2)^{3/2}}},

where G is Newton’s gravitational constant and x and y are the coordinates of the ball bearing in the plane perpendicular to the rod. The integral can be done in closed form and gives

F = {GMm\over{\sqrt{(x^2+y^2)(x^2+y^2+L^2/4)}}}.

Hence show that the equations of motion for the position x,y of the ball bearing in the xy-plane are

{d^2 x\over d t^2} = - GM {x\over r^2\sqrt{r^2+L^2/4}},
{d^2 y\over d t^2} = - GM {y\over r^2\sqrt{r^2+L^2/4}},

where r=\sqrt{x^2+y^2}.
Convert these two second-order equations into four first-order ones using the techniques of Section 8.3. Then, working in units where G=1, write a program to solve them for M=10, L=2, and initial conditions (x,y)=(1,0) with velocity of +1 in the y direction. Calculate the orbit from t=0 to t=10 and make a plot of it, meaning a plot of y against x. You should find that the ball bearing does not orbit in a circle or ellipse as a planet does, but has a precessing orbit, which arises because the attractive force is not a simple 1/r^2 force as it is for a planet orbiting the Sun.

8.13 Planetary orbits

This exercise asks you to calculate the orbits of two of the planets using the Bulirsch–Stoer method. The method gives results significantly more accurate than the Verlet method used to calculate the Earth’s orbit in Exercise 8.12.
The equations of motion for the position x,y of a planet in its orbital plane are the same as those for any orbiting body and are derived in Exercise 8.10 on page 361:

{d^2 x\over dt^2} = -GM {x\over r^3},
{d^2 y\over d t^2} = -GM {y\over r^3},

where G=6.6738\times10^{-11}\,\mathrm{m^3\,kg^{-1}\,s^{-2}} is Newton’s gravitational constant, M=1.9891\times10^{30}\,kg is the mass of the Sun, and r=\sqrt{x^2+y^2}.
Let us first solve these equations for the orbit of the Earth, duplicating the results of Exercise 8.12, though with greater accuracy. The Earth’s orbit is not perfectly circular, but rather slightly elliptical. When it is at its closest approach to the Sun, its perihelion, it is moving precisely tangentially (i.e., perpendicular to the line between itself and the Sun) and it has distance 1.4710\times10^{11}\,m from the Sun and linear velocity 3.0287\times10^4\,\mathrm{ms^{-1}}.

Write a program, or modify the one from Example 8.7, to calculate the orbit of the Earth using the Bulirsch–Stoer method to a positional accuracy of 1km per year. Divide the orbit into intervals of length H = 1 week and then calculate the solution for each interval using the combined modified midpoint/Richardson extrapolation method described in this section. Make a plot of the orbit, showing at least one complete revolution about the Sun.
Modify your program to calculate the orbit of the dwarf planet Pluto. The distance between the Sun and Pluto at perihelion is 4.4368\times10^{12}\,m and the linear velocity is
6.1218\times10^3\,\mathrm{ms^{-1}}. Choose a suitable value for H to make your calculation run in reasonable time, while once again giving a solution accurate to 1km per year.

You should find that the orbit of Pluto is significantly elliptical – much more so than the orbit of the Earth. Pluto is a Kuiper belt object, similar to a comet, and (unlike true planets) it’s typical
for such objects to have quite elliptical orbits.

Exercise 8.18 Oscillating chemical reactions

The Belousov–Zhabotinsky reaction is a chemical oscillator, a cocktail of chemicals which, when heated, undergoes a series of reactions that cause the chemical concentrations in the mixture to oscillate between two extremes. You can add an indicator dye to the reaction which changes color depending on the concentrations and watch the mixture switch back and forth between two different colors for as long as you go on heating the mixture.

Physicist Ilya Prigogine formulated a mathematical model of this type of chemical oscillator, which he called the “Brusselator” after his home town of Brussels. The equations for the Brusselator are

{dx \over dt} = 1 - (b+1)x + ax^2y
{dy \over dt} = bx - ax^2y.

Here x and y represent concentrations of chemicals and a and b are positive constants.

Write a program to solve these equations for the case a=1, b=3 with initial conditions x=y=0, to an accuracy of at least \delta=10^{-10} per unit time in both x and y, using the adaptive Bulirsch–Stoer method described in Section 8.5.6. Calculate a solution from t=0 to
t=20, initially using a single time interval of size H=20. Allow a maximum of n=8 modified midpoint steps in an interval before you divide in half and try again.

Make a plot of your solutions for x and y as a function of time, both on the same graph, and have your program add dots to the curves to show where the boundaries of the time intervals lie. You should find that the points are significantly closer together in parts of the solution where the variables are changing rapidly.

Hint: The simplest way to perform the calculation is to make use of recursion, as described in Exercise 8.17.