Flight of a projectile

• Download demo project - 8.0 Kb
• Download source - 4.3 Kb

Introduction

This article describes how to numerically solve a set of second order differential equations with initial conditions. Nyström modification of the fourth order Runge-Kutta method is explained first. Then the method is applied to two problems: to find the trajectory of a flying projectile and to calculate coupled oscillations of a mechanical system with two degrees of freedom. The article is intended as a continuation of my previous article about solitary differential equations.

Background

Many textbooks bring a simplified approach to the calculation of a flying projectile's path. The simplification is based on the assumption that the horizontal and vertical components of projectile motion are independent and can be treated separately. This enables the setup of two equations of motion: one for the vertical throw under the force of gravity and the other for the horizontal motion with no acceleration, which keeps that velocity constant. Each of these two equations can be solved individually yielding quite realistic results.

Observing the situation more closely, however, we find that some tie between the horizontal and vertical movements should exist. The reason is aerodynamic drag. For velocities higher than, say, 20 km/h, the drag of a moving body is proportional to the square of the velocity. If we suppose that the drag works directly against the velocity direction, then the square dependency causes coupling of the horizontal and vertical forces that act on the projectile. Imagine the ball moving with velocity v, like in the following picture.

Because of ball symmetry, the drag D is in the exactly opposite direction of velocity v. This means that the ratio of vertical and horizontal components of both velocity v and drag D should be the same.

v_y / v_x = D_y / D_x

Considering the square dependence of drag on velocity, we conclude that the vertical component of drag, D_y, relies both on the vertical and horizontal components of velocity v. The same is true for horizontal drag component. So, whenever we want to take the aerodynamic drag of a projectile into account we should solve a set of two tied equations of motion, not two individual equations. Let us now throw a ball with initial velocity v₀ and a given elevation angle.

At any following instant t, the horizontal and vertical forces acting on the ball are:

F_x = D_xF_y = G + D_y

where D_x is the horizontal and D_y the vertical component of aerodynamic drag. G is the ball weight. Obviously,

G = - m g

where m is the ball mass and g is the acceleration due the gravity. The air drag can be given as

D = 1/2 r_a v² c_d A

Here, r_a is the air density, c_d is the drag coefficient and A is the ball cross-sectional area. Total velocity is:

v = sqrt(v_x² + v_y²)

Drag components are then obtained as follows:

D_x = - D v_x / v
D_y = - D v_y / v

Dividing the resulting forces by the ball mass, we get two components of the ball's acceleration.

d²x/dt² = F_x / m
d²y/dt² = F_y / m

In other words, we obtained a set of two ordinal second order differential equations with initial conditions. There is a numerical method called Nyström modification of the fourth order Runge-Kutta method which suits well to such problems. The method claims that:

y_n+1 = y_n + h (dy/dt)_n + h²/6 (k₁ + k₂ + k₃)
(dy/dt)_n+1 = (dy/dt)_n + h/6 (k₁ + 2 k₂ + 2 k₃ + k₄)

where coefficients k₁ to k₄ are:

k₁ = f(t_n, y_n, (dy/dt)_n)
k₂ = f(t_n + 1/2 h, y_n + 1/2 h (dy/dt)_n + h²/8 k₁, (dy/dt)_n + h/2 k₁)
k₃ = f(t_n + 1/2 h, y_n + 1/2 h (dy/dt)_n + h²/8 k₂, (dy/dt)_n + 1/2 k₂)
k₄ = f(t_n + h, y_n + h (dy/dt)_n + h²/2 k₃, (dy/dt)_n + h k₃)

Here, f(t_n, y_n, (dy/dt)_n) is a function returning the value of the second derivative. It relies on present time t_n, function value y_n and the first derivative (dy/dt)_n. y_n+1 is the next function value, (dy/dt)_n+1 is the next first derivative and h is a time step. The nice feature of the Nyström method is that the same formulae can be applied both to one individual differential equation and also to a set of equations. In this second case, however, the function values y_n, first derivatives (dy/dt)_n, second derivatives (d²x/dt²)_n and even coefficients k₁ to k₄ are all vectors, not single numbers. This detail is the only challenge that we face moving from a numerical solution of one individual equation to a solution of a differential equation set. In the following section, we show how vectors can be represented in C# and how they keep our programming of the Nyström algorithm transparent.

Using the code

As we saw above, the algorithm to solve a set of second order differential equations works with vectors instead of single numeric variables. So, our first step should be to write a class suitable to hold and operate vectors.

PointN class

Because the word vector is used a bit informally among programmers, I call my container of n real numbers PointN. This name means "point in n-dimensional space." Here are member variables and two constructors:

public class PointN
{
    int n;
    double [] x;

    public PointN(int n)
    {
        this.n = n;
        x = new double[n];
    }

    public PointN(PointN a)
    {
        n = a.n;
        x = new double[n];
        for (int i = 0; i < n; i++)
        x[i] = a.x[i];
    }
}

The member variable n defines number of vector components and x is an array of real numbers. Each item of x holds one vector component. The first constructor simply sets the number of vector components and creates the array x. The second constructor is a copy constructor; it not only sets the number of components and creates x array, but also copies array values from the vector a. To get a quick approach to individual vector components, I define an indexer:

public double this[int i]
{
    get
    {
        if (i < 0 || i >= n)
            return 0;
        else
            return x[i]; 
    }
    set
    {
        if (i >= 0 && i < n)
            x[i] = value;
    }
}

Because my goal is to get simple and clear code for the Nyström integration algorithm, I define an addition of two vectors, a and b, and a multiplication of a vector a by a real number c:

public static PointN operator +(PointN a, PointN b)
{
    int n = a.n;
    PointN r = new PointN(n);
    for (int i = 0; i < n; i++)
        r.x[i] = a.x[i] + b.x[i];
    return r;
}

A subtraction can be done in the same manner. Now the multiplication of vector a by a real number c from the left:

public static PointN operator *(double c, PointN a)
{
    int n = a.n;
    PointN r = new PointN(n);
    for (int i = 0; i < n; i++)
        r.x[i] = c * a.x[i];
    return r;
}

The multiplication from the right side should also be added. Having readied our PointN class, we can now try to code an integration procedure.

Nyström integration

First, we should decide which shape the function calculating the second derivative should have. An interface is an ideal mean to do this:

public interface SecondDerivativeN
{
    PointN GetValue(double t, PointN y, PointN dy);
}

The interface is called SecondDerivativeN, which stands for "second derivative in n dimensions." It declares one method only, GetValue. The method takes time t, a vector of positions y and a vector of velocities dy and returns a vector of second derivatives. The integration is coded in NystromN class:

public class NystromN
{
    SecondDerivativeN secondDerivative;
    double t0;
    double h;
    PointN y0, dy0, d2y0;

    public NystromN(SecondDerivativeN secondDerivative, double t0,
        PointN y0, PointN dy0, double h, PointN d2y0)
    {
        this.secondDerivative = secondDerivative;
        this.t0 = t0;
        this.y0 = y0;
        this.dy0 = dy0;
        this.h = h;
        this.d2y0 = d2y0;
    }

    public int Step(out double t1, out PointN y1, 
        out PointN dy1, out PointN d2y1)
    {
        double h2 = h * h;

        PointN k1 = d2y0;
        PointN k2 = secondDerivative.GetValue(
            t0 + h/2, y0 + h/2 * dy0 + h2/8 * k1, dy0 + h/2 * k1);
        PointN k3 = secondDerivative.GetValue(
            t0 + h/2, y0 + h/2 * dy0 + h2/8 * k2, dy0 + h/2 * k2);
        PointN k4 = secondDerivative.GetValue(
            t0 + h, y0 + h * dy0 + h2/2 * k3, dy0 + h * k3);

        t1 = t0 + h;
        y1 = y0 + h * dy0 + h2/6 * (k1 + k2 + k3);
        dy1 = dy0 + h/6 * (k1 + 2 * k2 + 2 * k3 + k4);
        d2y1 = secondDerivative.GetValue(t1, y1, dy1);

        t0 = t1;
        y0 = y1;
        dy0 = dy1;
        d2y0 = d2y1;

        return 0;
    }
}

The constructor takes an instance of the SecondDerivativeN interface as its first argument. Other arguments are: initial time t0, initial vectors of positions and velocities y0 and dy0, time interval h and finally, an initial vector of second derivatives d2y0. This last one is a bit questionable. Why should we provide a value of the second derivative when the method calculating it was passed as the first argument? Well, this is my favorite way of organizing the integration process. The example of NystromN usage will illustrate this approach.

Projectile in flight

As an example, we use the problem mentioned in the Backgroud section. The wooden spherical ball having a rough surface (material density 600 kg/m³, diameter 0.1 m) is thrown under an elevation angle of 45 degrees. So, both the horizontal and vertical velocity components have the same value: 50 m/s in our case. Let the air density be 1.29 kg/m³ and the drag coefficient be 0.5.

The FlyingBall2 class takes all necessary parameters from its constructor and stores them in its member variables. dB is the ball diameter, roB the material density, cD the drag coefficient, g the acceleration due to gravity and roA the air density.

class FlyingBall2 : SecondDerivativeN
{
    double dB;
    double roB;
    double cD;
    double g;
    double roA;
    double mB;
    double aB;

    public FlyingBall2(double dB, double roB, double cD, double g, double roA)
    {
        this.dB = dB;
        this.roB = roB;
        this.cD = cD;
        this.g = g;
        this.roA = roA;

        double v = Math.PI * Math.Pow(dB, 3) / 6;
        mB = roB * v;
        aB = 0.25 * Math.PI * dB * dB;
    }

    PointN SecondDerivativeN.GetValue(double t, PointN y, PointN v)
    {
        double vx = v[0];
        double vy = v[1];
        double vxSqr = vx * vx;
        double vySqr = vy * vy;
        double vSqr = vxSqr + vySqr;
        double vTotal = Math.Sqrt(vSqr);

        double cosAlpha = 1;
        double cosBeta = 0;
        if (vTotal > 1.0e-12)
        {
            cosAlpha = vx / vTotal;
            cosBeta = vy / vTotal;
        }

        double drag = 0.5 * cD * roA * vSqr * aB;
        double dragX = -drag * cosAlpha;
        double dragY = -drag * cosBeta;

        PointN force = new PointN(2);
        force[0] = dragX;
        force[1] = -mB * g + dragY;

        PointN acceleration = force/mB;
        return acceleration;
    }
}

The SecondDerivateiveN.GetValue method implements the SecondDerivativeN interface. It calculates the two-dimensional vector of acceleration using the time t, the x and y coordinates in the y vector, and the v_x and v_y components of the velocity vector in v. The method determines the total velocity vTotal first. Then it gets two parameters that define the direction of velocity, cosAplha and cosBeta. Drag magnitude is stored in the drag variable and, with the help of ratios cosAplha and cosBeta, the drag components dragX and dragY are obtained. At this moment, we have all forces taking action on the projectile. Summing them with respect to their directions, we get the vector of the resulting force. The last step is to divide this vector by the projectile mass and we finally have the acceleration vector. By definition, acceleration is the second derivative of position, so the rest is now up to the integration routine NystromN.

static void Test1()
{
    FlyingBall2 projectileA = new FlyingBall2(0.1, 600, 0.0, 9.81, 1.29);
    SecondDerivativeN accelFunctionA = (SecondDerivativeN)projectileA;
    FlyingBall2 projectileB = new FlyingBall2(0.1, 600, 0.5, 9.81, 1.29);
    SecondDerivativeN accelFunctionB = (SecondDerivativeN)projectileB;

    double t0 = 0;
    double tmax = 12;
    double h = 0.5;

    PointN startPosition = new PointN(2);
    startPosition[0] = 0;
    startPosition[1] = 0;

    PointN startVelocity = new PointN(2);
    startVelocity[0] = 50;
    startVelocity[1] = 50;

    PointN startAccelerationA = accelFunctionA.GetValue(
        t0, startPosition, startVelocity);
    NystromN integratorA = new NystromN(
        accelFunctionA, t0, startPosition, startVelocity,
        h, startAccelerationA);
    PointN startAccelerationB = accelFunctionB.GetValue(
        t0, startPosition, startVelocity);
    NystromN integratorB = new NystromN(
        accelFunctionA, t0, startPosition, startVelocity,
        h, startAccelerationB);

    Console.WriteLine("TEST1 - Flight of projectile");
    Console.WriteLine(" t,s xA,m yA,m xB,m yB,m");
    Console.WriteLine("------------------------------------");
    Console.WriteLine("{0,4:F1}{1,8:F2}{2,8:F2}{3,8:F2}{4,8:F2}", t0,
        startPosition[0], startPosition[1], 
        startPosition[0], startPosition[1]);

    double t;
    PointN positionA, velocityA, accelerationA;
    PointN positionB, velocityB, accelerationB;
    do
    {
        integratorA.Step(out t, out positionA, 
            out velocityA, out accelerationA);
        integratorB.Step(out t, out positionB, 
            out velocityB, out accelerationB);
        Console.WriteLine("{0,4:F1}{1,8:F2}{2,8:F2}{3,8:F2}{4,8:F2}", t,
            positionA[0], positionA[1],
            positionB[0], positionB[1]);
    } while (t < tmax);
    Console.WriteLine();
}

To make things a bit more difficult, I decided to calculate two projectile trajectories: the first without air drag and the second having air drag included. Here, the results turned out to be really respectable. I declared two FlyingBall2 type objects: projectileA, which neglects the drag, and projectileB, which considers the drag. The difference is in the third argument of the constructor. In the case of projectileA, I passed the value of the drag coefficient as 0. For the projectileB object, I used a value 0.5, which is consistent with a sphere having a rough surface. The initial conditions given by the vectors startPosition and startVelocity were the same for both cases. Initial accelerations should differ; hence the need for vectors startAccelerationA and startAccelerationB. Certainly, vectors for calculated positions, velocities and accelerations should also be declared twice.

The integration is made in a do while cycle. Notice that the Step methods of both integrators are called in each cycle turn and that both of the calculated positions are displayed. Here are the results, first in table form:

Here are both trajectories graphically:

Coupled oscillations

Imagine the mechanical system with two degrees of freedom shown by the following picture.

The forces acting on masses m₀ and m₁ result from the deformations of springs k₀, k₁ and k₀₁. While the force produced by spring k₀ depends on the position of the m₀ mass only -- and the analogical rule is true of spring k₁ and mass m₁ -- the deformation of the spring k₀₁ is influenced by the positions of both masses. For this reason, the motions of both masses are interrelated. To describe this system, a set of equations of motion should be used. Here is a class that defines the accelerations of the two masses:

class CoupledOscillators : SecondDerivativeN
{
    double m0, m1;
    double k0, k1, k01;

    public CoupledOscillators(double m0, double m1, 
        double k0, double k1, double k01)
    {
        this.m0 = m0;
        this.m1 = m1;
        this.k0 = k0;
        this.k1 = k1;
        this.k01 = k01;
    }

    PointN SecondDerivativeN.GetValue(double t, PointN x, PointN v)
    {
        PointN force = new PointN(2);
        force[0] = -k0 * x[0] + k01 * (x[1] - x[0]);
        force[1] = -k1 * x[1] - k01 * (x[1] - x[0]);

        PointN acceleration = new PointN(2);
        acceleration[0] = force[0] / m0;
        acceleration[1] = force[1] / m1;
        return acceleration;
    } 
}

Member variables m0 and m1 hold masses m₀ and m₁, in kilograms. Variables k0, k1 and k01 store the stiffnesses of springs k₀, k₁ and k₀₁, suitable units being Newtons per meter. The method SecondDerivativeN.GetValue implements the SecondDerivativeN interface. The forces acting on both masses are calculated and stored in the force vector. Note the opposite signs in the contributions to shared spring k₀₁. This is the result of Newton's second law: If m₁ drags m₀ rightwards, then m₀ should push m₁ to the left. Eventually, the forces are divided by the masses and the two accelerations, which are in fact the second time-derivatives of the masses' positions. The integration is arranged similarly to our first example:

static void Test2()
{
    CoupledOscillators couple = new CoupledOscillators(
        1.0, 1.0,
        40.0, 40.0, 10.0);
    SecondDerivativeN accelFunction = (SecondDerivativeN)couple;

    double t0 = 0, tmax = 6, h = 0.05;
    PointN startPosition = new PointN(2);
    startPosition[0] = 0.1;
    PointN startVelocity = new PointN(2);
    PointN startAcceleration = accelFunction.GetValue(
        t0, startPosition, startVelocity);
    NystromN integrator = new NystromN(
        accelFunction, t0, startPosition, startVelocity,
        h, startAcceleration);

    Console.WriteLine("TEST2 - Coupled oscillation");
    Console.WriteLine(" t,s x0,m x1,m v0,m/s v1,m/s");
    Console.WriteLine("-------------------------------------------");
    Console.WriteLine("{0,4:F2}{1,10:F3}{2,10:F3}{3,10:F3}{4,10:F3}", t0,
        startPosition[0], startPosition[1],
        startVelocity[0], startVelocity[1]);

    double t;
    PointN position, velocity, acceleration;
    do
    {
        integrator.Step(out t, out position, out velocity, out acceleration);
        Console.WriteLine("{0,4:F2}{1,10:F3}{2,10:F3}{3,10:F3}{4,10:F3}", t,
            position[0], position[1],
            velocity[0], velocity[1]);
    } while (t < tmax);
    Console.WriteLine();
}

The member variable couple holds an object of CoupledOscillators type, which describes the behavior of the mechanical system. Parameters of the system are defined immediately in the couple constructor. startPosition and startVelocity define the initial conditions of the system and are passed to the integrator constructor, as well as the startAcceleration vector.

One step of the integration process is provided by the integrator.step method and is immediately followed by the print of results. A graph built from this print shows the exchange of energy between the two masses. At the beginning, the mass m₀ has its maximum amplitude and the mass m₁ is at rest. At the time of 2 seconds, the position of m₁ is at its maximum while m₀ moves slowly, close to its equilibrium location. As time passes, the cycle of energy exchange repeats.

Points of interest

For me, the most exciting feature demonstrated in this article is the fact that the C# code for the Nyström integration method can remain almost identical whether we're dealing with only one individual equation of motion or with a set of interrelated equations. Working with a set of equations, we used vectors in the role played by real numbers in the case of one standalone equation. To keep the new code as similar to the version with one equation as possible, we should define addition and multiplication operators for the vector type, which is not difficult.

Another remarkable item was the usage of two separate integrations in one turn of the cycle of the ball path calculation. One integration used a model without air drag while the other worked with equations that considered the drag. Thus, two different systems could be handled simultaneously, which simplified the output of results and their comparison.

History

• 8 June, 2007 -- Original version posted