физика — Реализация Рунге Кутты четвертого порядка в переполнении стека

Поэтому я должен начать с того, что это моя первая попытка выполнения исполняемого скрипта и мой первый скрипт на языке c ++. Имея это в виду, моя проблема заключается в следующем. Я пытаюсь создать сценарий, который объединит движение снаряда с использованием метода Рунге Кутты Четвертого Порядка. В настоящее время мой скрипт использует метод Эйлера и работает просто отлично. Однако, когда я пытаюсь реализовать свою функцию Runge Kutta, я получаю абсолютный мусор.

Например, с моей текущей интеграцией Эйлера я получаю следующее:

#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;
}

Затем я запускаю простой скрипт bash для тестирования с различными начальными углами

#!/bin/bash

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

Наконец, я отображаю свои результаты с помощью скрипта gnuplot.

#!/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'

Как уже упоминалось, моя интеграция Эйлера с кодом, показанным выше, приводит к следующему рисунку (который я знаю, чтобы быть правильным)

Тем не менее, когда я пытаюсь изменить свой код, чтобы использовать Runge Kutta, как я сказал абсолютный мусор.

#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;
}