Chapter Fifteen

Differential Equations with Java and Maple; Projectile Motion with Drag*

In this chapter we solve for projectile motion with drag by solving simultaneous, second-order ordinary differential equations (ODEs) using both Java and Maple. This is a fairly realistic problem that is often missing from undergraduate education, but whose solution looks familiar to students. The ODE solver used in our Java program is simple, useful for other problems as well, and provides some insight into numerical differentiation. The Maple solution, however, gets to be surprisingly involved at times due to the equations being both second-order and coupled, and having a solution that cannot be plotted simply. A comparison of the two approaches is quite educational. We have marked this chapter as optional, in part because there are no new materials that are used elsewhere, and in part because some students may not yet be familiar with differential equations. However the subject matter is important for many fields, and it does provide a good balance to the solution of a partial differential equation studied in Chapter 20.

15.1 PROBLEM: PROJECTILE MOTION WITH DRAG

We want to determine if the inclusion of air resistance in projectile motion leads to trajectories that are much steeper at their ends than their beginnings (golf balls that appear to drop out of the sky). We saw in Chapter 12, “Flow Control via Logic,” that frictionless projectile trajectories are parabolas, symmetric about their midpoints, and so do not describe balls dropping out of the sky.

In Figure 12.1 we showed trajectories of a projectile fired from a cannon without friction (solid curve) as well as one with air resistance or drag included (dashed curve). The projectile is fired at time t=0 with velocity V0 at an angle θ relative to the horizon. It remains in the air for a total or hang time t=T and hits the ground a horizontal distance or range x=R away from the origin, with all these quantities apparently changing if air resistance is included. The problem solution requires us to write a program that predicts or simulates these different trajectories for a given set of initial conditions.

15.2 MODEL: VELOCITY-DEPENDENT DRAG

The basic physics behind this problem is Newton’s second law of motion, which relates force, mass, and acceleration:

image

where the bold symbols indicate vector quantities. This is a second-order ordinary differential equation. A differential equation because it contains a derivative; an ordinary derivative because there is only one independent variable and so no partial derivatives (no ∂), the time t; and second-order because it is a second derivative. Because the equation of motion involves vectors, there are separate differential equations for the x and y components of each vector:

image

where we have divided through by the mass, switched the derivative to the LHS (a standard form), and substituted components for the frictional force Ff and the gravitational force.

The force of friction Ff is not a basic force of nature with a definite form for all situations. Indeed, it is just an approximate description of the physics of viscous flow, with no one expression being accurate for all velocities. We know it always opposes motion, which means it is in a direction opposite to that of the velocity. The simplest model, and the one often studied in texts [M&T 88], assumes that the force of air resistance is proportional to some power n of the projectile’s speed:

image

The v/|v| term ensures that the frictional force changes direction, to keep opposing motion, when the velocity changes sign. If n =1, or any odd power, then we may have the frictional force proportional to -vn and get the correct sign. However, if n is even, then we have to use these absolute values to ensure the correct sign. Though a frictional force proportional to a power of the velocity is a more accurate description than a constant force, it is still a simplification. Indeed, physical measurements indicate the power n appears to change with velocity, and so the most accurate model would be a numerical one that uses the empirical velocity dependence n(v).

With a simple power-law dependence, the equations of motion take the form:1

image

We shall consider three values for n, each of which represents a different model for the air resistance: (1) n = 1 for low velocities; (2) n = 1, for medium velocities; and (2) n =3/2 for high velocities. For comparison, we will compare these solutions with air resistance to the analytic results for the frictionless case:

image

15.3 ALGORITHM: NUMERICAL DIFFERENTIATION

In Chapter 5, “Solving Equations, Differentiation,” we used Maple to explore the accuracy of numerical differentiation. We now apply those techniques in the Java solution of differential equations. Most calculus courses start off with the definition of the derivative as:

image

It is a challenge to apply this definition on a computer, because its finite word length causes the numerator to fluctuate between 0 and round-off error as the denominator approaches zero. A simple and effective solution is to make h small but not zero:

Equation (15.8), known as the forward-difference rule, is not the most accurate approximation [CP 05], but it will suffice.

15.4 MATH: SOLVING DIFFERENTIAL EQUATIONS

When air resistance is present, the force on the projectile depends upon the value of the velocity, which means that the acceleration is not constant. This makes an analytic solution to the differential equation difficult, but it is not a challenge to a numerical solution. The trick is to figure out an approximation that takes the solution at the initial time 0 and advances it a small amount Δt. Because Δt is small, it is not hard to find an approximate solution. One then takes the solution so found and uses it as the new “initial conditions” to provide the solution at time 2Δt. The process is continued until the solution at the desired time is found. With [x(t = 0), y(t = 0)] and [vx(t = 0),vy(t = 0)] as the initial conditions, we have all the information we need to start the solution and will keep it going as long as we like.

The simplest algorithm for solving differential equations is a variation of Euler’s rule. It is based on the fact that for a sufficiently short time interval Δt, we may assume that the acceleration is constant and use the equations for constant acceleration from elementary physics to advance the velocity and position for a short time from t to t + Δt:

image

Here ax(t) is the variable acceleration at time t. Likewise, there are similar equations for y and vy. To make the algorithm simpler, we assume that Δt is so small that the terms quadratic in Δt are negligible with respect to the terms linear in Δt:

image

This means that as long as we keep the time step Δt small, we may ignore the acceleration term in (15.9). The acceleration still affects the solution via (15.10), but it does so by changing the velocity, which, in turn, changes the position. For this simple method to work, we must choose small values of Δt. As a rule of thumb, if the physical system under study has a characteristic time or period of T, then we would start with Δt ~ T/100 and then see the effect of decreasing it. A higher-order technique might require fewer steps but use more complicated formulas.

As our example, we take the simple n =1 law for air resistance:

image

The terms on the right-hand sides of these equations are the accelerations we need for the algorithm. The algorithm starts with the given initial position (the origin) and the velocities vx,o, vy,o, and then steps through many time steps. We label the time steps with the index i = 1…N. For each time step, the position and velocity will be updated according to the algorithm:

image

The acceleration due to gravity enters only for the y component of velocity, since gravity acts only in the y direction. And as we have said, the forces act by directly changing the velocity, but change the position only indirectly through the changes in velocity.

15.4.1 Implementation: ProjectileAir.java

The implementation of the algorithm is given in the code ProjectileAir.java shown in Listing 15.1. The class ProjectileAir contains:

1.  a main method that directs the computation;

2.  a plotNumeric method that solves the equations of motion numerically and plots the results by calling the PtPlot package; and

3.  a plotAnalytic method that calls the PtPlot package to plot the analytic solution for the frictionless case.

Examine how the PtPlot package is imported in line 2, before the class ProjectileAir is declared. On lines 4 and 5, after the class declaration but before any method declaration, the initial parameters for the problem are set as final variables, namely, constants that cannot be modified by any method. Because they are declared as static class variables, they will be accessible to all methods within the class, without having to give them as explicit arguments to the methods. However, methods in other classes do not have access to these static class variables.

Observe how on line 10 of ProjectileAir, we calculate the hang time T, height H, and range R for a frictionless projectile. These are useful to compare with the results from the simulation including air resistance. (If air resistance does not lead to a reduction in all three quantities, then we would suspect our simulation.) A plot object is created and labeled on lines 12–15. After that, data points are added to the plot object using the object-oriented notation myPlot.addPoint. The plot is drawn on your screen via the creation of a PlotApplication on line 19, and should look like the trajectory in Figure 12.1.

This program makes an indirect use of methods that is worthwhile noting. Normally a method is like a function that takes many arguments but returns only a single value for the function. Here we are repeatedly calling a method that keeps adding data points to an object without ever returning a value to the calling program. As discussed in Chapter 17, the exception to the one-value-returned rule occurs when an object or abstract data type (something with multiple parts) is used as an argument to a method. Because objects have multiple parts they are called by reference and not by value, which means that a reference to the location in memory where the object resides is passed to the method and not the value of different parts of the object. On account of this, the method may change the values of the different parts of the object (the plot in our case) as it resides in memory, as long as it does not change the location in memory where the object resides.

Observe next how the numerical algorithm used to solve the differential equation is contained in the method plotNumeric on line 21. It decides on the value of the time step dt by taking the total hang time for frictionless flight and dividing it by the input parameter N. It follows then that you make the time step dt small by making the number of time steps N large. We give some instructions on how to do that in the “Assessment” section. The rest of the method steps through time in small steps, uses Euler’s algorithm (15.10) to advance the position and velocity vector, and then adds points to the plot for each time step.

The plotAnalytic method on line 39 is similar to the plotNumeric method, in that it steps through time in small steps and add points to the plot for each time step. However, it computes the position and velocity vectors using the analytic expressions (15.5) for frictionless flight.

15.5 ASSESSMENT: BALLS FALLING OUT OF THE SKY?

1.  Compile and run ProjectileAir and check that you get a plot similar to that in Figure 12.1. If your program is not plotting, then you may need to review some of the materials on using PtPlot in Chapter 11, “Visualization with Java.”

2.  In general, it is not be possible to compare analytic and numerical results to realistic problems, because analytic expressions do not exist. However, we do know the analytic expressions (15.5) for the frictionless case, and we may turn friction off in our numerical algorithm and then compare the two. This is not a guarantee that we have handled friction correctly, yet it is a guarantee that we have something wrong if the comparison fails. Add some lines of code to ProjectileAir.java that calls the plotNumeric method a second time, this time with the friction coefficient k=0.

3.  We have deliberately given you a program in which the total number of time steps N is rather small, and, inversely, the time step dt is rather large. Consequently, the results from the given program will be of low precision. Repeat the calculation with ever-increasing values for N, until you cannot tell the difference on the plot between the numeric and the analytic results for frictionless flight.

4.  Once the results look good to you, determine the level of precision of the calculation by comparing the actual numbers for the analytic and the numerical calculations. If necessary, increase N until the relative difference is less than 1%. Use this value of N for future computations.

5.  Add lines to the program so that it prints out and plots up the components of the velocity as a function of time.

6.  Study your plots of the velocity components versus time and draw a line on them that appears to be the terminal velocity for each.

7.  Print out several plots showing the effect of air resistance for different initial angles θ. For each, determine the value of R, T, and H, and compare the three with the corresponding values for frictionless flight. All three should decrease.

8.  We know that for frictionless flight, the maximum range R occurs for θ= 45°, and that there is the same range for trajectories with θ or 90° — θ (for example, for 30° and 60°).

9.  Investigate how both these statements change when drag is included. In other words, if you want to hit your golf ball with a maximum range, should you hit it at an angle that is greater than or less than 45°? 9. Investigate and describe how the range R varies with initial velocity vQ for a fixed value of θ. Explain the results you obtain (does it go further, for example?).

10.  Up until this point you have looked at results for the model of friction (15.3) with n =1. This corresponds to a low-velocity model. Modify Projectile.java to investigate n = 2, (high-velocity model) and n = 3/2 (medium-velocity model). To make a realistic comparison, adjust the value of k for the latter two cases such that the initial force of friction, kVn, is the same for all three cases.

Listing 15.1 ProjectileAir.java

image

15.5.1 Improved Algorithm: Verlet*

In developing the Euler integration algorithm (15.10), we ignored the direct affect of the acceleration on the position, as is present in (15.9). We now want to modify the algorithm so that it includes the acceleration term directly in the solution for x(t + Δt). An approach that is nearly as simple as Euler yet appears to work very well is called the Verlet algorithm. It uses second-order terms for xn, but only first-order ones for vn:

image

We see that while the algorithm for the velocity is only first-order in Δt, it makes up for this somewhat by using the average acceleration over the time interval and not just its value at the beginning of the interval. The only problem with this method is that it is not fully self-contained; to determine vn+1 you first have to use the Euler method to obtain an+.

Exercise: Implement the Verlet algorithm and compare the number of steps needed to obtain the three-place precision with both the Euler and Verlet methods. (The better algorithm should use fewer steps.) ♠

15.6 MAPLE: DIFFERENTIAL-EQUATION TOOLS

Maple is competent at finding solutions to most any differential equation that you give it. It many cases it can find an analytic solution with no help from you, while in other cases you may have to give it some hints based on your understanding of the mathematics of your equation. If all else fails, or if all you want is a plot of the solution, ask Maple to solve the differential equation numerically.

In this section we solve the projectile-with-drag problem in Maple. Because we have already done this with Java, it affords you the opportunity to compare the two approaches to the same problem. You should decide for yourself; we find the Java solution easier and faster.

Warning: Although the commands we have given here appear to work properly on various versions of Maple, we have noted that the order in which Maple returns the solution changes with different versions. To cite an instance, [y(t), x(t)] in one case and [x(t), y(t)] in another. Consequently, you may have to reverse the order of some of the arguments in the plotting commands given below.

We start by repeating the differential equations for projectile motion in a uniform gravitational field with no friction, and their solution:

image

where we have assumed x(0) = y(0) = 0, and used the symbols vo, x and vo, y to denote the x and y components of the initial velocity v0. We start by solving these equations with Maple and seeing if we get these same solutions. To start simply, first we solve for just the y motion, and then for the simultaneous x and y motions. Then we add in air resistance. This careful approach to a final solution is standard in scientific problem solving.

To keep from repeatedly entering complicated expressions, we give our equations names and enter the names rather than long expressions. We write our differential equations using the diff command, and later show how to do it with the D operator (see Chapter 5 for our discussion of derivatives). Yet even if diff is used for the equations, D must be used for the initial conditions.

There are various ways to enter differential equations, and we demonstrate several of them (all of which give the same answer of course). We use the variables diff eq and diff eq2 for two different representations of the equation for the y component of a projectile’s motion, and solve each with the dsolve command:

image

Two simpler ways of entering the same differential equation is with the D operator:

image

We solve this equation with the dsolve. Even though the command accepts a number of arguments, some of which direct Maple along specified paths for solutions, we leave it all up to Maple:

image

We see that we do indeed get the general solution with integration constants C1 and C2 to be determined by the initial conditions (the leading underscores indicate that the computer has chosen these names).

We incorporate the initial conditions into the solution by setting C1 = V y0 and C2 = 0, or by giving dsolve the initial conditions along with the equations and letting it do all the work for us:

image

This looks great and is the right answer, yet we have had to be careful. First, even if the differential equation is entered with diff, the initial conditions must be entered with the D. Second, because we have a second-order ODE, there are two initial conditions (displacement and velocity), and we group them as the set init_con, with items separated by a comma, but order irrelevant.

As a check of the solution, we take the derivative of y(t) to obtain an expression for the velocity, and then take the derivative of the velocity to see if we get the acceleration back. To do this, we need to keep in mind that Maple returns an expression in the form of an equation for y(t) as the solution. Even though the expression looks like a function, it is not. Even though it is possible to transform y(t) into a function, the procedure is not straightforward. (One needs to first define a procedure that contains the solution, and then create a function that calls the procedure). We take a simpler approach in which we define the expression soltn to contain the solution, and then manipulate soltn:

image

image

15.6.1 Extract Right-Hand Side: rhs

In order to plot the solution, all parameters must be assigned numerical values. In addition, since the expression for the solution is an equation with y(t) on the LHS, we must use Maple’s rhs() function to extract the RHS. We extract the velocity and acceleration from the solutions by taking time derivatives:

image

Check over how the velocity starts off positive and then gets more and more negative, and that the acceleration is constant at g = 9.8 in the negative y direction.

15.6.2 Second-Order ODE with Plot

Maple has the DEplot command that both solves a differential equation and plots up the solutions. As before, while our solution was analytic and contained initial conditions as parameters, to plot a graph we need to have a completely numeric answer, and so we must assign values to all parameters. The value of the parameter g is assigned first, and the values of the initial conditions are placed within the DEplot command:

image

image

Observe how we include the initial conditions as a list for the third argument.

15.6.3 System of ODEs

When the systems in an environment depend on each other, we usually have simultaneous differential equations to solve. As a case in point, in our projectile problem we have simultaneous equations for the x and y motions (we solved only for the y motion so far). The same commands and procedures are used for a set of equations as for a single equation. The additional step is to define the variable that gave the equation to solve to now be a set of equations (“set” because the order does not matter):

image

Here we use D to define the equations, although diff would do. The solution is obtained by specifying the set of equations {diff_eqs} as the first argument to dsolve and the set,{x(t), y(t)} the solution we desire, as the second argument:

image

We see that Maple does find the expected forms of the solution (15.16), with the constants still to be determined. As before, we include the initial conditions into the set of equations we give as the first argument to dsolve, this time with four conditions:

image

Indeed, we get the familiar solutions as a set within braces. If we want numeric values for x and y, we must assign numeric values to the parameters.

15.7 MAPLE SOLUTION: DRAG ∝ VELOCITY

We now return to our projectile-with-drag problem, applying some of the Maple tools we have just learned. If the frictional force is proportional to the first power of the velocity, the equations to solve are:

image

Because the x and y motions are independent, there is no mixing of the x motion into the y equations, and vice versa, and we could actually solve each equation separately. However, simultaneous equations generally are coupled, and we will solve them as if they were. We enter the equations into dsolve as a set with elements separated by commas. We use a set because order does not matter, and later we will place brackets around the elements. We define a variable diff_eqs as the set of equations, and then solve the set. Because the equations involve the second derivatives, there will be a diff with respect to t$2 and t:

image

Observe that we must indicate the t dependence of x and y as x(t) and y(t), because otherwise Maple would treat them as constants whose derivatives vanish.

Now we see if Maple is smart enough to find an analytic solution. We enter the set of equations as the first argument to dsolve with braces around diff_eqs to indicate that it is a set:

image

We see that Maple returns a set of solutions with four constants. Good, we are on our way. Rather than solving for the constants in terms of the initial conditions, we define a set of initial conditions, and include them as a second argument to dsolve. We specify initial velocities as the initial values for the first derivatives, and use the D() operator (diff does not work for initial conditions):

image

image

This appears to be the analytic solution we want, and it is worth investigating its properties. Yet before we do that, we need to see if the velocities exhibit reasonable behaviors, that is, if they attain a terminal speed and then stop increasing. We do that by taking the derivative of our solution:

image

Indeed, we see that the x component of velocity decreases exponentially as a function of time from its initial value and approaches zero for long times. This is a consequence of the frictional force always opposing the velocity in the x direction, but with no other force present to increase the velocity in the x direction. In contrast, the y velocity approaches a terminal value -g/k for long times. This is a consequence of the force of gravity F = -mg exactly balancing the frictional force F = mkv when v = -g/k:

image

Once the parameters are assigned, we check the solution:

image

We see our solution in its numeric form, as a set with a comma separating the equation for x(t) from that of y(t). Though they make the equations hard to read, the zeros after the decimal point indicate the Maple’s level of precision for this floating-point calculation (it became floating point when we made g and k floatingpoint numbers).

15.8 EXTRACT OPERANDS

To plot up the results we need to extract the RHSs of both equations. For a simple expression this was done with the rhs command. Now that soltn is a set of equations, we need to indicate which equation it is in the set that we want the RHS of. For this purpose Maple has the op command to extract operands (pieces) from an expression. Consequently, rhs(x) is equivalent to op(2, x). First we test that we know how to extract each equation, and then we take RHSs:

image

At last we have what we need to make our plots. We plot the position and velocity (time derivative) by entering lists (in square brackets) as the first argument to the plot command (a list because order matters for the plot):

image

These plots show an initial linear increase of x(t) with time t, as expected for constant velocity. However, as time increases, the action of the frictional force in the x direction is to keep decreasing Vx until, for large time, Vx vanishes (except that the ground gets in the way). Then the x motion ceases and x remains constant. In turn, y(t) shows a parabolic dependence on time initially, as expected for uniform acceleration, but only a linear increase for later times. The linear increase of y with time indicates that vy has attained its terminal value as the drag force exactly balances that of gravity. As expected from the analytic expressions, and our discussion in the previous paragraph, we see that Vx approaches zero and Vy approaches a constant value for large times (or would if the ground did not get in the way).

The usual plot of a trajectory is a plot of y(t) as the ordinate versus x(t) as the abscissa. As you will recall from our discussion of visualization in Chapter 4, this type of plot is known as a parametric plot. A parametric plot is constructed by moving the time limits into the list that is used as the first argument (the list contains {x(t), y(t)} and the range of t):

image

We see a trajectory that is much distorted from the symmetric parabola that occurs for frictionless flight, and looks similar to that computed with Java and shown in Figure 12.1. In fact, this looks very much like the type of trajectory seen for a golf ball or a baseball in which the ball rises at first and then appears to “drop out of the sky” at the end of its trajectory.

15.8.1 Solution for R and T

We now want to determine if we can solve for the range R and hang time T for the general projectile-with-friction problem. We start by again placing our solution in symbolic form:

image

The range R corresponds to the value of x at which the height y =0 (in addition to the initial firing). We start our solution for the hang time T at which y(T) = 0 by extracting x(t) and y(t):

image

This just tells us the initial time, which we know, but is not what we want. Instead, let us be more specific and ask for the time at which the height y just becomes negative, namely, when it just hits the ground:

> solve(Y<0,t>);

We see that Maple returns nothing. It apparently has given up. If we look at the expression for y(t), we see that it contains the time t in both an exponential and a linear term. Apparently, setting Y =0 leads to a transcendental equation that has no analytic solution. The only solution, then, is obtained numerically, and we shall do that with Java (although Maple will also work).

15.9 DRAG ∝V2 (EXERCISE)

Modify the analysis performed for a drag force proportional to the first power of the velocity so that it is appropriate for a force proportional to the square of the velocity. Maple should be able to solve this analytically as well, and you should be able to follow all of the steps we have for Model 1. Note: since v2 does not change sign when v does, you will need to add a v/|v| factor in the frictional force for the y motion (it is not needed for the x motion because the x component of velocity is always positive).

15.10 DRAG ∝V3/2

We now want to solve the projectile-with-friction problem for a drag force proportional to the 3/2 power of the velocity. The equations we need to solve are:

image

The acute reader might notice that there is a problem with this last equation; namely, if vy is negative, then the square root operation, which is part of raising vy to the 3/2 power, will return an imaginary number. This is clearly not what we want. We solve that problem by taking the absolute value of vy before raising it to a power:

image

This expression is simpler to read and easier to enter into Maple. (An alternative approach would be to use the sign function, which extracts the sign of its argument, to keep track of the sign of the velocity and adjust the frictional force appropriately.)

We set about solving this model just as we did the problem with drag proportional to the first power. First we define the set of differential equations and give it the name diff eqs:

image

This looks fine, and so we go on to look for a general solution by first entering the initial conditions:

image

Warning: You may have to stop the following command by hand. We now ask Maple to solve the differential equations labeled diff_eqs with the initial conditions labeled InitConds:

image

We see from the time it takes Maple to search for a solution, from the volume of output that Maple produces, and from the RootOf command being returned, that Maple is having trouble finding a simple solution to our problem. In any case, an analytic solution with this level of complexity is not illuminating, and probably not even good for computation (the numerous subtractions and evaluations of the multivalued tan-1 and ln functions is error-prone).

Now we try a numerical approach. We start it by assigning values to the initial conditions and parameters:

image

We get a numeric solution using the same procedures as before, only now by including the type=numeric argument to dsolve:

image

Regardless of this not appearing to be a clear signal that a solution is in hand, this is Maple’s way of telling us that it has written a procedure named soltn that will produce a numerical solution to the equations. Recall, a procedure in Maple is like a function that requires more than one line to define. It is like a method in Java or a subroutine in Fortran. In the present case the argument to the procedure tells us that the procedure takes an argument and then returns the solution. Since Maple assumes we are solving for y(x), the argument is indicated generically as x. In our case, the second argument to the solve command was the list [x(t), y(t)], and so we are solving for x and y as functions of the time t. Therefore the argument to the procedure is the time t. (The subscript rkf45 to x is meant for aficionados; it indicates that a fourth order Runga-Kutta numerical method is used, and that it uses adaptive step size to obtain fifth-order precision.)

We now have a procedure soltn(t) that numerically solves the differential equations for a single input value of time t and returns [x(t), y(t)] in some form. To see how this works, let us look at the solution returned for time 0 (which should just be the initial conditions we entered):

image

We use Maple to pick out the individual pieces of the solution:

image

image

15.10.1 Defining Functions from Procedures

It seems like we should now be able to plot up the solution. However, Maple’s plot command accepts expressions and functions as input but not procedures. (We discussed Maple programming and procedures in Chapter 8.) This means we cannot have plot call our procedure and plot it. It is not much work to define a function from a procedure and then to plot up the function; essentially, you just define an arrow function as a call to your procedure. Here we do that for the x and y positions and velocities:

image

We see similar results to Model 1, only now with an even more drastic “drop out of the sky effect” at the end of the trajectory.

image

Figure 15.1 Left: The gravitational force on a planet a distance r from the sun. The x and y components of the force are indicated. Right: Output from the applet PlanetRHL showing the precession of a planet’s orbit when the gravitational force 1/r4.

15.11 EXPLORATION: PLANETARY MOTION*

Newton’s explanation of the motion of the planets in terms of a universal law of gravitation is one of the great achievements of science. He was able to prove that the planets traveled along elliptical paths with the sun at one vertex, and with periods that agree with observation. All Newton needed to postulate is that the force between a planet of mass m and the sun of mass M is given by

image

where r is the distance between the planet of mass m and sun of mass M, and G is a universal constant. The minus sign indicates that the force is always attractive and lies along the line connecting the planet and sun, as indicated on the left of Figure 15.1. The hard part for Newton was solving the resulting differential equations, since he had to invent calculus to do it. Whereas the analytic solution is complicated, the numerical solution is not.

Even for planets, the basic equation of motion is

image

with the force now given by (15.20). If we look at Figure 15.1 we see that

image

with F given by (15.20). If we write the equation of motion in component form and substitute (15.22) for the force components, we obtain the differential equations we need to solve:

image

15.11.1 Implementation: Planet.java

On the CD you will find the applet Planet. We suggest that you run the class file to see how the solution behaves. To keep the calculation simple without changing the physics, assume that the units we use are such that GM = 1, and that the initial conditions are [Feyn 63]:

image

1.  Modify projectileAir.java so that it solves (15.23) as functions of time.

2.  Make the number of time steps large enough so that the planet’s orbit repeats on top of itself.

3.  You may need to make the time step small enough so that the orbit closes upon itself (it should) and just repeats. This should be a nice ellipse.

4.  Experiment with initial conditions until you obtain the ones that produce a circular orbit (a circle is a special case of an ellipse).

5.  Once you have good precision, see the effect of progressively increasing the initial velocity until the orbit opens up and becomes an hyperbola.

6.  Use the same initial conditions as produced the ellipse. This time investigate the effect of the power in (15.20) being 1/r4 rather than 1/r2. You should find that the orbital ellipse now rotates (precesses), as in Figure 15.1.

15.12 KEY WORDS

equations of motion Euler’s method forward difference

numerical ODE solution numerical differentiation ODE

15.13 SUPPLEMENTARY EXERCISES

1.  Use Maple, Java, or both to find the solution of the differential equation

image

with the initial condition n(0) = 100. Plot up the answer for λ = 0.3.

2. In Chapter 16 you will encounter an RLC electrical circuit and the differential equation describing it. We consider here the same circuit without a capacitor and described by the differential equation

image

Solve this equation for the current in the circuit for the same values of R and L as in the problem, and for a constant V (t) (same magnitude as in problem). Plot your solution.

3.  Consider exponential growth starting with an initial population N(0) = 2.

a. Solve the equation describing exponential growth

image

Assign the integration constants to correspond to n(0) = 100, and λ = 0.3, and plot up the solution. b. Solve for exponential growth with a modulated growth parameter

image

This is the differential-equation version of the logistics map studied in Chapter 19, “Discrete Math, Arrays as Bins.” The term in square brackets is insignificant for small N(t).

4. Solve for and plot up the solution to the differential equation

image

with initial condition y(0) = 1.

a.  Verify that Maple cannot find an analytic solution to this equation.

b.  Have Maple find a numerical solution to this equation.

5.  Solve for and plot up the solution to the differential equation

image

subject to the initial conditions y(0) = 1, y'(0) = 0. Hint: If Maple cannot find an analytic solution, you may need to try a numerical one.

1  For more realistic models used for friction, these equations may become coupled, e.g., because (V 2)x = (Vx)2. For simplicity, we ignore the coupling of the equations, although the technique for the numerical solution remains unchanged.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset