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 RungeKutta 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_{0} and a given elevation angle.
At any following instant t, the horizontal and vertical forces acting on the ball are:
F_{x} = D_{x}F_{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^{2} c_{d} A
Here, r_{a} is the air density, c_{d} is the drag coefficient and A is the ball crosssectional area. Total velocity is:
v = sqrt(v_{x}^{2} + v_{y}^{2})
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^{2}x/dt^{2} = F_{x} / m
d^{2}y/dt^{2} = 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 RungeKutta method which suits well to such problems. The method claims that:
y_{n+1} = y_{n} + h (dy/dt)_{n} + h^{2}/6 (k_{1} + k_{2} + k_{3})
(dy/dt)_{n+1} = (dy/dt)_{n }+ h/6 (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})
where coefficients k_{1} to k_{4} are:
k_{1} = f(t_{n}, y_{n}, (dy/dt)_{n})
k_{2} = f(t_{n} + 1/2 h, y_{n} + 1/2 h (dy/dt)_{n} + h^{2}/8 k_{1}, (dy/dt)_{n} + h/2 k_{1})
k_{3} = f(t_{n} + 1/2 h, y_{n} + 1/2 h (dy/dt)_{n} + h^{2}/8 k_{2}, (dy/dt)_{n} + 1/2 k_{2})
k_{4} = f(t_{n} + h, y_{n} + h (dy/dt)_{n} + h^{2}/2 k_{3}, (dy/dt)_{n} + h k_{3})
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^{2}x/dt^{2})_{n} and even coefficients k_{1} to k_{4} 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 ndimensional 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^{3}, 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^{3} 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.0e12)
{
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 twodimensional 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:
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 m_{0} and m_{1} result from the deformations of springs k_{0}, k_{1} and k_{01}. While the force produced by spring k_{0} depends on the position of the m_{0} mass only  and the analogical rule is true of spring k_{1} and mass m_{1 } the deformation of the spring k_{01} 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_{0} and m_{1}, in kilograms. Variables k0
, k1
and k01
store the stiffnesses of springs k_{0}, k_{1} and k_{01}, 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_{01}. This is the result of Newton's second law: If m_{1 }drags m_{0} rightwards, then m_{0} should push m_{1} to the left. Eventually, the forces are divided by the masses and the two accelerations, which are in fact the second timederivatives 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_{0} has its maximum amplitude and the mass m_{1} is at rest. At the time of 2 seconds, the position of m_{1} is at its maximum while m_{0} 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