-4
Итак, я должен начать с того, что это моя первая попытка исполняемого скрипта и мой первый скрипт C++. Учитывая это, моя проблема заключается в следующем. Я пытаюсь создать скрипт, который должен интегрировать движение снаряда с использованием метода Рунге Кутты четвертого порядка. В настоящее время мой скрипт использует метод Эйлера и работает отлично. Однако, когда я пытаюсь реализовать свою функцию Runge Kutta, я получаю абсолютный мусор.Реализация четвертого ордера Runge Kutta в C++
Например, с моей текущей интеграцией Эйлера я получаю следующее:
#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;
}
тогда я бег через простой Баш скрипт для тестирования с различными начальными углами
#!/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'
Как уже упоминалось, интеграция Эйлера с t код, показанный выше, дает следующий рисунок (который, как я знаю, является правильным)
Однако, когда я пытаюсь изменить свой код, чтобы использовать 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;
}
Перед публикацией прочитайте [* «Как спросить» *) (http://stackoverflow.com/help/how-to-ask). –