Here's a dip into one of the more complex branches of
mathematics, Analysis - the area of mathematics that deals largely with
change and rates of change. I'll warn readers in advance that this is not
my area of expertise, so what is presented here is my simple-minded
understanding of some background information.

#### Introduction

Differential equations contain information about the rate of
change of an dependent variable with respect to some independent variable.
To solve the equation, we integrate (sum) these little changes to get a final
answer. Sometimes there are analytical solutions to these types of
problems, often there are not. In these cases computers come to the
rescue with numerical solutions.

#### Runge-Kutta

If the equation meets certain conditions then the
classical 4th order Runge-Kutta method is a good choice for computer solution
because it is both fast and accurate. The "4th order" part
refers to the fact that the algorithm takes a weighted average of 4
estimates of the derivative for each calculated point which reduces total error
in proportion to the 4th power of the chosen step size. One of the
conditions for Runge-Kutta is "**Ordinary**", only one independent
variable, and the abbreviation **ODE** is frequently used as shorthand
reference for **Ordinary Differential Equation**. We'll do that from
now on. Many real world phenomena can be described by ODEs
that meet the conditions required for Runge-Kutta including our current
interest, time/location problems. I'll assume that we're in the
time/location domain for the rest of this discussion.

#### Initial Conditions

ODEs do not completely specify where an object will be at some
point in time unless we know where it was and how fast it was moving at some
other point in time. Thus we also need to specify **Initial Conditions**
when using Runge-Kutta. (A variation for future implementation
uses Boundary Values instead - we know where it object was at two points in
time, but not it's velocity.)

#### 2nd Order

The usual notation for derivatives uses** '** marks, so for
variable **x**, **x'** is the first derivative and **x''** is the
second derivative.

Runge-Kutta requires that ODEs be **linear**, that is contain
first derivatives only. Fortunately it can handle
systems with multiple equations and multiple dependent variables and it is easy
to split an equation that contains a 2nd derivative into two equations that
contain only first derivatives by assigning a new variable (and equation) **x2=x'**.
Now any prior occurrence of **x''** becomes **x2'.** In
our second order cases, the called procedure will take care of estimating the
first derivative given a second derivative, as well as estimating the variable
value from the 1st derivative. It will be up to us to
tell it what 2nd derivative value to use at any given point in time, T,
given **T**, **X**, and **X'** values.

#### Systems of Equations

Sometimes there are "coupled" systems where several
objects are affect each other, planets in the solar system, or springs that
connect things together for example. Described below is a second
Runge-Kutta procedure that can solve these systems if we provide an array of
initial conditions, one for each element on the system, and an array of
functions to do the necessary 2nd derivative calculations at any point in
time.

Enough background, let's move on to -

### The Implementation

Unit **URungeKutta4** defines a type **Float** (currently
double) that applies to all of the floating point variables used.

There are two procedures presented here:

**RungeKuttaIC**

**RungeKuttaIC** solves a single *Second Order ODE
with Initial Conditions*. The calling sequence is

RungeKutta2ndOrderIC( LowerLimit,
UpperLimit, InitialValue, InitialDeriv,
ReturnInterval, CalcInterval, Error, UserFunc, CallBackFunc);

where

Lower Limit and Upper
Limit are type Float and give the range of
the independent variable (0 to 2.5 seconds for example)

InitialValue and** **InitialDeriv**
**are the Float values of the independent
variable when the dependent variable is at LowerLimit.

ReturnInterval is the Float
independent variable specifying the independent variable increment
for result points to be passed back to the
caller.

CalcInterval is the
Float independent variable increment between calculation points
Note: ReturnInterval should be a multiple of CalcInterval.
CalcInterval influences the accuracy of the result so we may, for example
calculate results every hundredth of a second (CalcInterval) but only return
result to the caller for each tenth of a second (ResultInterval).

Error:
return error codes

0 - All OK

1 - ReturnInterval must be greater than 0.

2 - ReturnInterval must be greater than CalcInterval.

3 - UpperLimit cannot be the same as Lowerlimit.

UserFunc is the name of a caller
supplied function with calling sequence : Function
UserFunc(T, X, XPrime:Float):Float; which must return a second
derivative value for the (T, X, Xprime) values received. T is the current
independent variable value, X the current dependent variable value, and XPrime
the current 1st derivative value for X. This function is defined with the
"**of object" **condition,** **so** **your function should be
defined within an object, usually your **TForm1** main form object.

CallBackFunc is the name of a
caller supplied function which receives the T, X, and XPrime values every
ReturnInterval increments of the independent variable. Calling
sequence is Function CallBackFunc(T, X, XPrime:Float):
boolean; Parameters are the same as described above - return **True**
to continue the calculation, return **False** to stop the Runge-Kutta
procedure before the independent variable reaches Upperlimit.

#### RungeKutta2ndOrderIC_System

Solves a __system__ of *Second Order ODEs with
Initial Conditions*. The calling sequence is

RungeKutta2ndOrderIC_System( LowerLimit,
UpperLimit, InitialValues, ReturnInterval, CalcInterval,
Error, NumEquations, UserFuncs, Test2CallBackFunc);

where

Lower Limit and Upper
Limit are type Float and give the range of
the independent variable (0 to 2.5 seconds for example)

InitialValues is defined as type TnVector;
an array of [1..10] of TnData records containing
the initial conditions for each dependent variable**. **TnData**
**defined as :

TNData = record

X : Float;

xPrime : Float;

end;

ReturnInterval is the Float
independent variable increment between result points to be passed back to the
caller.

CalcInterval is the
Float independent variable increment between calculation points
Note: ReturnInterval should be a multiple of CalcInterval.

Error:
return error codes

0 - All OK

1 - ReturnInterval must be greater than 0.

2 - ReturnInterval must be greater than CalcInterval.

3 - UpperLimit cannot be the same as Lowerlimit.

NumEquations is the number of
equations in the system.

UserFuncs is type TFuncVect,
an array [1..10] of functions to be called while calculating points The calling
function for each function is : Function
UserFuncN(v:TNVector):Float; which must return a second derivative
value for the Nth dependent variable. V is an
array [0..10] of TnData records as described above.
Function type is defined with the "**of object" **condition**, **so**
**your functions should be defined within an object, usually your **TForm1**
main form object. Note that variables are numbered 1 through 10; V[0].X
is used to contain the current value of the independent variable.

CallBackFunc is the name of a
caller supplied function which receives an array of Tndata records containing
the T, X, and XPrime
values each dependent variable every ReturnInterval
increments of the independent variable. Calling sequence is Function
Test2CallBackFunc(v:TNVector):boolean; Parameters are the same as
described above - return **True** to continue the calculation, return **False**
to stop the Runge-Kutta procedure before the independent variable reaches Upperlimit.

### Running/Exploring the Program

Whew! Examining sample code will probably be an easier way
to learn how to use these routines. The code below contains two test cases from
Borland's Numerical Methods Toolbox, a simple "spring-mass" system and
"two pendulums coupled by a weak spring". Also a Simple Pendulum
I coded myself just to make sure understood the ODE part. I
added some simple animation "just for fun".