12,559,868 members (47,977 online)
alternative version

49.6K views
51 bookmarked
Posted

# Flight of a projectile

, 8 Jun 2007
 Rate this:
Numerical solution of a set of second order differential equations

## 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.

vy / vx = Dy / Dx

Considering the square dependence of drag on velocity, we conclude that the vertical component of drag, Dy, 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 v0 and a given elevation angle.

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

Fx = DxFy = G + Dy

where Dx is the horizontal and Dy 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 ra v2 cd A

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

v = sqrt(vx2 + vy2)

Drag components are then obtained as follows:

Dx = - D vx / v
Dy = - D vy / v

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

d2x/dt2 = Fx / m
d2y/dt2 = Fy / 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:

yn+1 = yn + h (dy/dt)n + h2/6 (k1 + k2 + k3)
(dy/dt)n+1 = (dy/dt)n + h/6 (k1 + 2 k2 + 2 k3 + k4)

where coefficients k1 to k4 are:

k1 = f(tn, yn, (dy/dt)n)
k2 = f(tn + 1/2 h, yn + 1/2 h (dy/dt)n + h2/8 k1, (dy/dt)n + h/2 k1)
k3 = f(tn + 1/2 h, yn + 1/2 h (dy/dt)n + h2/8 k2, (dy/dt)n + 1/2 k2)
k4 = f(tn + h, yn + h (dy/dt)n + h2/2 k3, (dy/dt)n + h k3)

Here, f(tn, yn, (dy/dt)n) is a function returning the value of the second derivative. It relies on present time tn, function value yn and the first derivative (dy/dt)n. yn+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 yn, first derivatives (dy/dt)n, second derivatives (d2x/dt2)n and even coefficients k1 to k4 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/m3, 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/m3 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 vx and vy 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:

 Time, s Drag neglected Drag included Distance, m Height, m Distance, m Height, m 0.0 0.00 0.00 0.00 0.00 0.5 25.00 23.77 23.81 22.59 1.0 50.00 45.10 47.62 42.72 1.5 75.00 63.96 71.44 60.40 2.0 100.00 80.38 95.25 75.63 2.5 125.00 94.34 119.06 88.41 3.0 150.00 105.86 142.87 98.73 3.5 175.00 114.91 166.69 106.60 4.0 200.00 121.52 190.50 112.02 4.5 225.00 125.67 214.31 114.98 5.0 250.00 127.38 238.12 115.50 5.5 275.00 126.62 261.94 113.56 6.0 300.00 123.42 285.75 109.17 6.5 325.00 117.76 309.56 102.32 7.0 350.00 109.66 333.37 93.03 7.5 375.00 99.09 357.18 81.28 8.0 400.00 86.08 381.00 67.08 8.5 425.00 70.61 404.81 50.42 9.0 450.00 52.70 428.62 31.32 9.5 475.00 32.32 452.43 9.76 10.0 500.00 9.50 476.25 -14.25 10.5 525.00 -15.78 500.06 -40.72 11.0 550.00 -43.51 523.87 -69.63 11.5 575.00 -73.69 547.68 -101.00 12.0 600.00 -106.32 571.49 -134.83

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 m0 and m1 result from the deformations of springs k0, k1 and k01. While the force produced by spring k0 depends on the position of the m0 mass only -- and the analogical rule is true of spring k1 and mass m1 -- the deformation of the spring k01 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 m0 and m1, in kilograms. Variables `k0`, `k1` and `k01` store the stiffnesses of springs k0, k1 and k01, 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 k01. This is the result of Newton's second law: If m1 drags m0 rightwards, then m0 should push m1 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 m0 has its maximum amplitude and the mass m1 is at rest. At the time of 2 seconds, the position of m1 is at its maximum while m0 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

A list of licenses authors might use can be found here

## Share

 Software Developer Czech Republic
Software developer, mainly in Visual C++

## You may also be interested in...

 View All Threads First Prev Next
 Good start. JimOnDotNet12-Jun-07 5:28 JimOnDotNet 12-Jun-07 5:28
 Re: Good start. cokerms12-Jun-07 10:13 cokerms 12-Jun-07 10:13
 Re: Good start. JimOnDotNet12-Jun-07 12:17 JimOnDotNet 12-Jun-07 12:17
 Re: Good start. jussac13-Jun-07 2:20 jussac 13-Jun-07 2:20
 Last Visit: 31-Dec-99 18:00     Last Update: 27-Oct-16 19:37 Refresh 1