// RK++NET.cpp : main project file.
#include "stdafx.h"
using namespace System;
#include <map>
#include <vector>
#include <iostream>
const double gravity = -9.81;
template<typename T>
inline
T spring(const T& pos)
{
return -100.0 * pos;
}
template<typename T>
inline
T drag(double k, const T& pos)
{
return -k * pos;
}
template<typename A>
class EulerImplicit
{
public:
EulerImplicit(const A& accel):
accel(accel)
{
}
template<typename T>
inline
void operator()(T& pos, T& speed, const double delta) const
{
speed += delta * accel(pos, speed);
pos += delta * speed;
}
private:
const A& accel;
};
template<typename A>
class RK4
{
public:
RK4(const A& accel):
accel(accel)
{
}
template<typename T>
inline
void operator()(T& pos, T& speed, const double delta) const
{
const double half_delta = 0.5 * delta;
std::pair<T, T> k1 = accel(pos, speed);
std::pair<T, T> k2 = accel(pos + half_delta * k1.first, speed + half_delta * k1.second);
std::pair<T, T> k3 = accel(pos + half_delta * k2.first, speed + half_delta * k2.second);
std::pair<T, T> k4 = accel(pos + delta * k3.first, speed + delta * k3.second);
pos += delta / 6.0 * (k1.first + 2.0 * (k2.first + k3.first) + k4.first);
speed += delta / 6.0 * (k1.second + 2.0 * (k2.second + k3.second) + k4.second);
}
private:
const A& accel;
};
template<typename T>
class Force
{
public:
Force() : m_drag(1.0)
{}
T operator()(const T& pos, const T& speed) const
{
return drag(m_drag, speed) + spring(pos) + gravity;
}
private:
const double m_drag;
};
template<typename A>
class Adapt
{
public:
Adapt(const A& accel):
accel(accel)
{
}
template<typename T>
inline std::pair<T, T> operator()(const T& pos, const T& speed) const
{
std::pair<T, T> ret;
ret.first = speed;
ret.second = accel(pos, speed);
return ret;
}
private:
const A& accel;
};
template<typename T, typename I>
inline
void simulate(const I& intg_func, T& pos, T& speed, double t0, const double t_fin, const double delta)
{
while(t0 < t_fin)
{
t0 += delta;
intg_func(pos, speed, delta);
}
}
typedef std::vector< std::pair<double, double> > Tuples;
int main(array<System::String ^> ^args)
{
const double t0 = 0.0;
const double t_fin = 1000.0;
const double delta = 0.025;
Tuples initial_states;
initial_states.reserve(1000);
for (int i=0; i<1000; ++i)
initial_states.push_back( std::make_pair((double)i, 0.0) );
Tuples states(initial_states);
System::DateTime start0 = System::DateTime::Now;
Force<double> accel;
EulerImplicit< Force<double> > euler(accel);
for (Tuples::iterator it = states.begin();
it != states.end();
++it)
{
simulate(euler, it->first, it->second, t0, t_fin, 0.25 * delta);
}
System::TimeSpan diff0 = System::DateTime::Now - start0;
std::cout << "Time: " << diff0.TotalSeconds;
for (int i = 0; i < 20; ++i)
{
std::cout << " state = (" << states.at(i).first
<< "; " << states.at(i).second << ")" << std::endl;
}
states = initial_states;
start0 = System::DateTime::Now;
Adapt< Force<double> > adapted(accel);
RK4< Adapt< Force<double> > > rk4(adapted);
for (Tuples::iterator it = states.begin();
it != states.end();
++it)
{
simulate(rk4, it->first, it->second, t0, t_fin, delta);
}
diff0 = System::DateTime::Now - start0;
std::cout << "Time: " << diff0.TotalSeconds;
for (int i = 0; i < 20; ++i)
{
std::cout << " state = (" << states.at(i).first
<< "; " << states.at(i).second << ")" << std::endl;
}
Console::ReadKey(true);
return 0;
}
