实现Runge Kutta第四顺序在c ++ [英] Implementation of Runge Kutta Fourth Order in c++
本文介绍了实现Runge Kutta第四顺序在c ++的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
所以我应该开始说这是我第一次尝试一个可执行脚本，和我的第一个c ++脚本。承担这在我的心中我的问题如下。我试图创建一个脚本，它将使用四阶Runge Kutta方法集成弹丸运动。目前我的脚本使用的是Euler方法，工作正常。但是当我尝试实现我的Runge Kutta函数时，我得到绝对的垃圾。
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.
例如对于我当前的Euler集成，我得到以下：
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;
}
然后我运行一个简单的bash脚本测试各种初始角度
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
最后，结果到gnuplot脚本
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'
$ b b
如上所述，我的Euler与上面显示的代码集成生成下面的图（我知道是正确的）
As mentioned my Euler integration with the code shown above produces the following figure (which I know to be correct)
但是，当我尝试更改我的代码使用Runge Kutta就像我说的绝对垃圾。
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;
}
推荐答案
数字odeint，请参见此包含示例的网页以及这些标题
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
和例如这个Runge Kutta第4阶总览页面。
这篇关于实现Runge Kutta第四顺序在c ++的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
