Runge-Kutta algorithm as a C/C++ function: Passing functions
to functions
P. Gorham, updated 2/27/2018
1. Introduction
Because of the Runge-Kutta calculation, this week's lab is
kind of tedious and tricky to program. It is of great advantage if
we can generalize the form of the algorithm and create a separate function
that will do the job of the RK calculation for us. Then we can just
call the Runge-Kutta function to evolve our differential equation
solution, by passing to the RK function just the force laws it needs
to complete the calculation at each step. Then if we want to use it
again with a different force law, we will be good to go, with a
general purpose tool, at least for first-order differential
equations.
Assume that the variables xt and vt
will indicate the "updated" values of the position x(t) and velocity v(t) based on the RK changes,
and xtold, vtold are the original values x(t0) and v(t0) before we update
them. We would like our code to look something like this
pseudocode:
for(t=t0; t<Tmax;t += dt){
xt = xtold + RK(xarg1,xarg2,xarg3...); // the
new position is old one + RK increment
vt = vtold + RK(varg1,varg2,varg3...); //
the new velocity is old one + RK increment
xtold = xt; // update the old position
before going on to next one
vtold = vt; // and update the old
velocity also
cout << .....; // print out the
current position and velocity at each step
}
where the RK( ) function calculates the x and v increments according
to the general prescription for the Runge-Kutta methods (either 2nd
or 4th order) that we found in the lecture and lab this week.
Since we may want to do this with other force laws than the harmonic
oscillator, as we noted above, we want to have our RK( )
function independent of the type
of force law that we use. (Remember: the force law is the
right-hand-side (RHS) of the dv/dt equation in our problem.)
So, ideally, we would like one of the arguments to RK( ) to contain
the form of the function which calculates the f(t,y), which is the general
form for the RHS force law. We therefore pass the function an
argument which is another function -- the force law itself.
This can be done quite easily, but we first have to know how to
pass functions as arguments to another function. This
is a useful skill to develop in our C++ programming, so we
will take some time this week to understand it.
Here is an external
link which also discusses this topic, it may be helpful to
get another view of it.
2. General form
The general form for a function definition in which one of the
arguments is another function is easiest to understand in an
example. In the case of the Runge-Kutta method, we will be trying to
create a pair of functions that have the general form like this:
double RK(int ytype, double (*func1)( double , double) ,
double (*func2)( double , double), double
arg2 , double arg3, ....)
{
double tmp, tmp1, tmp2, retval
if(ytype==0) tmp = arg2 * arg3 *
(*func1)(tmp1,tmp2) // test if x or v
if(ytype==1) tmp = arg2 * arg3 *
(*func2)(tmp1,tmp2)
return (retval);
}
double func1( double darg1, double darg2)
{
return( darg1*darg2);
}
double func2( double darg1, double darg2)
{
return( darg1*darg2);
}
The function RK( ) takes
as input two pointers
to other functions called func1(
) or func2( ) and
these functions themselves take in a pair of doubles and return
another double. The function RK(
) also takes other double arguments, and it multiplies them
by the value returned by func1( ),
or func2 ( ) and
then itself returns another double. Whether we use func1( ) or func2( ) depends on an
external flag ytype, which indicates which of the canonical
variables, x or v we are calling the RK
function to get an update for.
Note that, although the usage of the passed function in RK( ) is a bit different than
you are used to, the declaration of the functions passed to to
it, func1( ) &
func2( ), are not affected by this different usage.
Also, the different usage applies ONLY to the case where the
function is passed as an argument. If you were to use the function func( ) itself in another
context where it was not passed as an argument to another function,
you would use it the usual way, by calling it with its own normal
arguments to get its normal return value.
3. A simple example: numerical derivative of a user-supplied
function
To illustrate how a function is passed to another function, we use
the following fairly simple program, deriv.cpp
which calculates a numerical derivative (or slope) of a user-defined
function at any given point x
and differential interval dx.
The program is not a very useful one, since the function must be
entered and recompiled each time it is changed, but it does illustrate the
method for passing functions as arguments. Take a look at it and
read through the comments.
// deriv.cpp:
// take the scalar numerical derivative of
// a user defined function
// P. Gorham March 2013 for Physics 305
using namespace std;
#include <iostream>
#include <iomanip>
#include <fstream>
#define USE_MATH_DEFINES
#include <cmath>
#include <cstdlib>
//---------------------FUNCTION DEFINITIONS HERE --------------------
double userfunc(double x) // this is the function whose derivative we want
{
double f_of_x;
// put your function in here
// use quadratic equation for example
f_of_x = x*x+4.*x-1.; // Right hand side is user supplied
return(f_of_x);
}
//----------------------------------------
/* this function of a function takes the numerical slope at point xval, for interval dx */
double derivative( double (*userfn)(double), double xval, double dx )
{
double dy_dx;
// next line gives derivative at position x according to:
// dy/dx(x) = (y(x+dx)-y(dx))/dx
dy_dx = ((*userfn)(xval+dx) - (*userfn)(xval))/(2.*dx);
return(dy_dx);
}
//-----------------MAIN HERE ---------------------------------
main(int argc, char *argv[])
{
double x,dx,dy_dx;
if(argc<2){
cerr << "usage; deriv [x value] [dx value]"<<endl;
exit(0);
}
x = atof(argv[1]);
dx = atof(argv[2]);
dy_dx = derivative(userfunc,x,dx);
cout << "At x= << x <<" dx= "<< dx << " derivative = "<< dy_dx << endl;
}
The function userfunc( )
is passed as a pointer to
a double-precision function. This is the reason for the
funny notation in the definition of the derivative( ) function below. Since derivative( ) is meant to
accept any function of one double-precision argument, the name of
its input function is NOT the same as the actual function--it is
just a dummy argument, which here is given the name userfn( ). The (*userfn) syntax
ensures that the address
of the actual function to be evaluated is given to derivative( ) as its
argument. This in turn means that any function can be placed in that
same address without changing the way that derivative( ) works.
You should try to understand this program completely before looking
at how this method applies to the Runge-Kutta calculation in the
next section.
4. Example programs for Runge-Kutta function arguments
An example of how this is implemented for the 2nd order Runge Kutta
case of a 1-dimensional force law (the harmonic oscillator or simple
pendulum, for example) is given in FRK2xv.cpp
.
/* FRK2xv: a function to calculate second-order Runge-Kutta
for 1-dimensional F = dp/dt = f(x(t),v(t)) (Newton force law)
the user must supply the functions f_x, f_v which determine
the return value.
Both x and v increments must be calculated at each step in
the integration of the equation, and xold and vold must be
reset after each step
--P. Gorham: updated 3/19/2007 for Physics 305, add explicit time dependence,
correct error in previous version */
double FRK2xv(int ytype, double (*f_x)(double, double, double),
double (*f_v)(double, double, double),
double t,double xold,double vold,double dt)
/* ytype = 0 for x(t) increment, 1 for v(t) increment,
xold, vold, previous position and velocity
t=time variable, dt = step size
double (*f_x)(double, double ,double) ==> pointer to a function of
three doubles that returns a double. Note the order: the f_x function
comes first in the argument list.
*/
{
double k1x, k1v, k2;
k1x = dt*(*f_x)(t, xold, vold);
k1v = dt*(*f_v)(t, xold, vold);
if(ytype==0){
k2 = dt*(*f_x)( t, xold+k1x/2.0, vold+k1v/2.0 );
}else{ /* ytype = 1 */
k2 = dt*(*f_v)( t, xold+k1x/2.0, vold+k1v/2.0 );
}
return(k2); // add this value to the previous v or x
}
/* the following are examples of the form of the functions needed.
In practice you should put the relevant functions in the source
code that contains the loop for integration of the position
and velocity */
/*the dx/dt function only requires returning the velocity--no action
double f_x (double x, double v)
{
return (v);
}
// Here is the actual force law
double f_dvdt( double x, double v)
{
insert the proper force law into the return value
double retval;
return( retval );
}
*/
Linking to object files: Using
your makefile more effectively
To use this in another program you either need to import all of the
source code into that program, or you can use the object file
FRK2xv.o and link it into your other program, which then allows you
to call the functions contained within FRK2xv.o.
To do this you will need to have the following lines in your
makefile, assuming you want to link FRK2xv.o with a program
harmoscRKr1.cpp (which is available
also for download as an example of how to use this technique):
FRK2xv.o:
g++ -c
FRK2xv.cpp
harmoscRKr1: FRK2xv.o
g++ $(GFLAGS)
harmoscRKr1.cpp FRK2xv.o -o harmoscRKr1 $(LIBS)
This construction tells the makefile interpreter that your
executable program harmoscRKr1 depends on the existence of FRK2xv.o
in the same folder as harmoscRKr1.c, and if FRK2xv.o is not present,
it will be created according to the instructions in the FRK2xv.o:
line.
Here is the program harmoscRKr1.cpp:
/* simple harmonic oscillator using external Runge-Kutta function */
/* P. Gorham, originally 3/6/2003 for UH Physics 305, updated 3/14/2007*/
/* updated 3/17/07 to use modified Runge-Kutta code that repaired
a subtle bug: FRK2xv.c --PG */
// -- updated to c++ 3/10/2013 --PG
using namespace std;
#include <iostream>
#include <iomanip>
#include <fstream>
#define USE_MATH_DEFINES
#include <cmath>
#include <cstdlib>
#define Tmax 50 // seconds
extern double FRK2xv(int, double (*)(double, double, double),
double (*)(double, double, double),
double ,double ,double ,double );
double f_x(double, double, double), f_v(double, double, double);
main(int argc, char *argv[])
{
double xt0,t0,vt0,vt,xt,t,dt,xtold,vtold;
ofstream outfile;
outfile.open("oscRK.dat");
xt0 = 1.0; // initial position
t0 = 0.0; // initial time
vt0 = 0.0; // initial velocity
dt = 0.01;
xtold = xt0; // set & print the initial conditions
vtold = vt0;
t = t0;
outfile << t << " " << xt0 << " " << vt0 << endl;
/* Here is the loop that propagates the motion:
vt is the new time, vtold the previous step time;
xt is the new position, xtold the previous;
after each step is calculated, the old is set
to the new, and the cycle is repeated */
for(t=t0; t<Tmax; t+= dt){
xt = xtold + FRK2xv(0,f_x,f_v,t,xtold,vtold,dt);
vt = vtold + FRK2xv(1,f_x,f_v,t,xtold,vtold,dt);
xtold = xt;
vtold = vt;
outfile << t << " " << xt << " " << vt << endl;
}
}
/* the following is just the formal definition: f(x) = dx/dt == v */
double f_x(double t, double x, double v)
{
return(v);
}
/* this is where the action is: should contain the force law: dv/dt = f(t,x,v) */
double f_v(double t, double x, double v)
{
double k = 1.0, m=1.0;
return( -k/m*x );
}
Notice the new type of declaration before main( ) in this program,
which uses the "extern" statement since the
FRK2xv function is now in an external module which will be linked at
the time you compile.
Try this program and modify the previous plot file to see how things
have improved. Try out some different
intervals (the "dt" variable) and see how large dt can be before
things start to get noticeably worse.
BACK to Runge Kutta Method