Implementation of Runge Kutta Fourth Order in c++

Question

So I should start by saying that this is my first attempt a executable script, and my first ever c++ script. Bearing this in my mind my problem is as follows. I'm attempting to create a script which shall integrate projectile motion with the use of the Fourth Order Runge Kutta method. Currently my script is using the Euler method and is working just fine. However when I attempt to implement my Runge Kutta function I get absolute garbage out.

For example with my current Euler integration I get the following:

#include <iostream> // Statndard Input Output
#include <stdlib.h> // Standard Library atof
#include <cmath> // Use Of Math Functions
#include <fstream> // File Stream Input Output
#include <string> // String Manipulation c_str
#include <sstream> // Used For Variable String Name

/* Allows for abbrevaited names
** without having to use std::cout
** for example we may simply type
** cout
*/
using namespace std; 

// Required Function Delcarations
double toRadians(double angle);
double xVelocity(double vel,double angle);
double yVelocity(double vel,double angle);
double rc4(double initState, double (*eqn)(double,double),double now,double dt);
double updatePosX(double currentPos,double vel,double deltaT);
double updatePosY(double currentPos,double vel,double deltaT);
double updateVelY(double yVel,double deltaT);
double updateVelX(double xVel,double deltaT);

int main(int argc, char *argv[]) //In Brackets Allows Command Line Input
{
    /* atof Parses the C string str, interpreting its 
    ** content as a floating point number and 
    ** returns its value as a double. */
    double v0 = atof(argv[1]);  
    double theta = atof(argv[2]);
    double timeStep = atof(argv[3]);
    double duration = atof(argv[4]);

    /* Simple printed message to the 
    ** console informing the user
    ** what set of initial conditions
    ** are currently being executed
    */
    cout << "Running Simulation" << endl;
    cout << "Velocity: " << v0 << " m/s" << endl;
    cout << "Angle: " << theta << " deg" << endl;
    cout << endl;

    //Setting x and y velocity 
    double vx = xVelocity(v0,toRadians(theta));
    double vy = yVelocity(v0,toRadians(theta));

    //Initial Conditions 
    double xPos = 0;
    double yPos = 0;
    double time = 0;

    /*Creating File Name With Vars
    ** Note: stringsteam is a stream
    ** for operating on strings. In 
    ** order to concatinate strings
    ** or produce statements like 
    ** those if java, we must use
    ** this stream. We CANNOT 
    ** simply concationate strings
    ** and variables
    */
    stringstream ss; 
    ss << "v0_" << v0 << "_ang_" << theta << ".txt";
    string fileName = ss.str();

    //Opening File For Writing
    ofstream myFile;
    myFile.open(fileName.c_str()); //Must be c style string

    // Checking Status of Attempt To Open   
    if(myFile.fail())
    {
        cout << "Failed To Open File For Writing" << endl;
    }

    else
    {   
        //myFile << "x pos \t y pos \t vx \t vy" << endl;

        // Doing Required Integration
        while (time <= duration && yPos >=0)
        {

            vx = updateVelX(vx,timeStep);
            vy = updateVelY(vy,timeStep);

            xPos = updatePosX(xPos,vx,timeStep);
            yPos = updatePosY(yPos,vy,timeStep);
            cout << "x Pos: " << xPos <<"\t y Pos: " << yPos << endl;

            time+=timeStep;

            myFile << xPos << "  " << yPos << "  " << vx << "  " << vy << endl; 

        }
    }

    //Closing File After Finished   
    myFile.close();

    return 0;
}

/* This function shall take a given
** angle in degrees and convert it
** to radians
*/
double toRadians(double angle)
{
    return (M_PI * (90-angle)) / 180.0;
}

/* This function shall take the inital
** angle at which the projectile is 
** launched and return the x componet
** of its velocity 
*/
double xVelocity(double vel,double angle)
{
    return vel * sin(angle);
}

/* This function shall take the inital
** angle at which the projectile is 
** launched and return the y componet
** of its velocity 
*/
double yVelocity(double vel,double angle)
{
    return vel * cos(angle);
}

/* This function shall be
** the X position of our 
** projectile
*/
double updatePosX(double currentPos,double vel,double deltaT)
{
    return currentPos + vel*deltaT;
}

/* This function shall be
** the Y posistion of our
** projecticle
*/
double updatePosY(double currentPos,double vel,double deltaT)
{
    return currentPos + vel*deltaT;
}

/* This function shall update
** the Y component of our 
** projectile's velocity 
*/
double updateVelY(double yVel,double deltaT)
{
    double g = 9.81;
    return yVel - g*deltaT;
}

/* This function shall update
** the X component of our
** projecticle's velocity
*/
double updateVelX(double xVel,double deltaT)
{
    return xVel;
}

/* This function shall be the fourth
** order Runge Kutta integration method
** and shall take as input y0 and return
** y+1.
**
** initState: Inital state of function
** (*eqn): Equation to be integrated
** now: Initial time to start integration
** dt: Timestep
*/
double rc4(double initState, double (*eqn)(double,double),double now,double dt)
{
    double k1 = dt * eqn(initState,now);
    double k2 = dt * eqn(initState + k1 / 2.0, now + dt / 2.0);
    double k3 = dt * eqn(initState + k2 / 2.0, now + dt / 2.0);
    double k4 = dt * eqn(initState + k3, now + dt);

    return  initState + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
}

I then run through a simple bash script to test with various initial angles

#!/bin/bash

for i in `seq 30 5 60`
do  
    ./projectile 700 $i 0.1 500
done

Lastly I plot my results to with a gnuplot script

#!/usr/bin/gnuplot

set terminal pdf color solid
set output "test.pdf"
set yrange[0:20000]
set xrange[0:55000]
plot for [i=30 : 60 : 5] 'v0_700_ang_'.i.'.txt' using 1:2 with lines title  'Angle '.i.' deg' 

As mentioned my Euler integration with the code shown above produces the following figure (which I know to be correct)

However, when I attempt to change my code to use Runge Kutta like I said absolute garbage.

#include <iostream> // Statndard Input Output
#include <stdlib.h> // Standard Library atof
#include <cmath> // Use Of Math Functions
#include <fstream> // File Stream Input Output
#include <string> // String Manipulation c_str
#include <sstream> // Used For Variable String Name

/* Allows for abbrevaited names
** without having to use std::cout
** for example we may simply type
** cout
*/
using namespace std; 

// Required Function Delcarations
double toRadians(double angle);
double xVelocity(double vel,double angle);
double yVelocity(double vel,double angle);
double rc4(double initState, double (*eqn)(double,double),double now,double dt);
double updatePosX(double currentPos,double vel,double deltaT);
double updatePosY(double currentPos,double vel,double deltaT);
double updateVelY(double yVel,double deltaT);
double updateVelX(double xVel,double deltaT);

int main(int argc, char *argv[]) //In Brackets Allows Command Line Input
{
    /* atof Parses the C string str, interpreting its 
    ** content as a floating point number and 
    ** returns its value as a double. */
    double v0 = atof(argv[1]);  
    double theta = atof(argv[2]);
    double timeStep = atof(argv[3]);
    double duration = atof(argv[4]);

    /* Simple printed message to the 
    ** console informing the user
    ** what set of initial conditions
    ** are currently being executed
    */
    cout << "Running Simulation" << endl;
    cout << "Velocity: " << v0 << " m/s" << endl;
    cout << "Angle: " << theta << " deg" << endl;
    cout << endl;

    //Setting x and y velocity 
    double vx = xVelocity(v0,toRadians(theta));
    double vy = yVelocity(v0,toRadians(theta));

    //Initial Conditions 
    double xPos = 0;
    double yPos = 0;
    double time = 0;

    /*Creating File Name With Vars
    ** Note: stringsteam is a stream
    ** for operating on strings. In 
    ** order to concatinate strings
    ** or produce statements like 
    ** those if java, we must use
    ** this stream. We CANNOT 
    ** simply concationate strings
    ** and variables
    */
    stringstream ss; 
    ss << "v0_" << v0 << "_ang_" << theta << ".txt";
    string fileName = ss.str();

    //Opening File For Writing
    ofstream myFile;
    myFile.open(fileName.c_str()); //Must be c style string

    // Checking Status of Attempt To Open   
    if(myFile.fail())
    {
        cout << "Failed To Open File For Writing" << endl;
    }

    else
    {   
        //myFile << "x pos \t y pos \t vx \t vy" << endl;

        // Doing Required Integration
        while (time <= duration && yPos >=0)
        {

            vx = rc4(vx,updateVelX,time,timeStep);
            vy = rc4(vy,updateVelY,time,timeStep);

            xPos = updatePosX(xPos,vx,timeStep);
            yPos = updatePosY(yPos,vy,timeStep);
            cout << "x Pos: " << xPos <<"\t y Pos: " << yPos << endl;

            time+=timeStep;

            myFile << xPos << "  " << yPos << "  " << vx << "  " << vy << endl; 

        }
    }

    //Closing File After Finished   
    myFile.close();

    return 0;
}

/* This function shall take a given
** angle in degrees and convert it
** to radians
*/
double toRadians(double angle)
{
    return (M_PI * (90-angle)) / 180.0;
}

/* This function shall take the inital
** angle at which the projectile is 
** launched and return the x componet
** of its velocity 
*/
double xVelocity(double vel,double angle)
{
    return vel * sin(angle);
}

/* This function shall take the inital
** angle at which the projectile is 
** launched and return the y componet
** of its velocity 
*/
double yVelocity(double vel,double angle)
{
    return vel * cos(angle);
}

/* This function shall be
** the X position of our 
** projectile
*/
double updatePosX(double currentPos,double vel,double deltaT)
{
    return currentPos + vel*deltaT;
}

/* This function shall be
** the Y posistion of our
** projecticle
*/
double updatePosY(double currentPos,double vel,double deltaT)
{
    return currentPos + vel*deltaT;
}

/* This function shall update
** the Y component of our 
** projectile's velocity 
*/
double updateVelY(double yVel,double deltaT)
{
    double g = 9.81;
    return yVel - g*deltaT;
}

/* This function shall update
** the X component of our
** projecticle's velocity
*/
double updateVelX(double xVel,double deltaT)
{
    return xVel;
}

/* This function shall be the fourth
** order Runge Kutta integration method
** and shall take as input y0 and return
** y+1.
**
** initState: Inital state of function
** (*eqn): Equation to be integrated
** now: Initial time to start integration
** dt: Timestep
*/
double rc4(double initState, double (*eqn)(double,double),double now,double dt)
{
    double k1 = dt * eqn(initState,now);
    double k2 = dt * eqn(initState + k1 / 2.0, now + dt / 2.0);
    double k3 = dt * eqn(initState + k2 / 2.0, now + dt / 2.0);
    double k4 = dt * eqn(initState + k3, now + dt);

    return  initState + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
}

Answer 1

These should be part of Boost Numeric odeint, see this pages with examples as well as these headers

boost/numeric/odeint/stepper/runge_kutta4.hpp
boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp
boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp
boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp
boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp

and eg this Runge Kutta 4th order overview page.