Runge-Kutta O.D.E. Algorithm

[Home]   [Puzzles & Projects]    [Delphi Techniques]   [Math topics]   [Library]   [Utilities]

 

Search

Search WWW

Search DelphiForFun.org

As of October, 2016, Embarcadero is offering a free release of Delphi (Delphi 10.1 Berlin Starter Edition ).     There are a few restrictions, but it is a welcome step toward making more programmers aware of the joys of Delphi.  They do say "Offer may be withdrawn at any time", so don't delay if you want to check it out.  Please use the feedback link to let me know if the link stops working.

 

Support DFF - Shop

 If you shop at Amazon anyway,  consider using this link. We receive a few cents from each purchase.   Thanks.


Support DFF - Donate

 If you benefit from the website,  in terms of knowledge, entertainment value, or something otherwise useful, consider making a donation via PayPal  to help defray the costs.  (No PayPal account necessary to donate via credit card.)  Transaction is secure.

Contact

Feedback:  Send an e-mail with your comments about this program (or anything else).

Search DelphiForFun.org only

 

 

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 xx' 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 variableTnData 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".  

 

Original Date: August 28, 2002

Modified:February 18, 2016

 
  [Feedback]   [Newsletters (subscribe/view)] [About me]
Copyright 2000-2017, Gary Darby    All rights reserved.