Download: https://github.com/dimitrs/odes
Below, listing 1. is an implementation of a generic algorithm – the Runga-kutta method which is well known for its use in the approximation of solutions of ordinary differential equations. Although the Runga-kutta method is well known and relatively simple, the generic algorithm implementation below is difficult to read. I won't attempt to explain how it works; I wrote it several months ago and I don't remember all the details. The reason why I list it here is to be able to compare it with a non-generic implementation of the algorithm I wrote in listing 2.
Listing 2. is far easier to read and understand. It is almost a carbon copy of the mathematical equations for the Runga-kutta4 method (see listing 3). The non-generic implementation is easy to read (assuming you are familiar with the equations) and it was very easy to write. I was done in about 30 minutes, whereas the generic version took me a couple of days to complete. So why make the effort ? The generic version allows you to use any container that supports forward iterators e.g. std::array, std::vector, C array, boost::ublas vectors, Blitz++. With the non-generic version you are tied into using only boost::ublas vectors.
As far as performance is concerned, the generic version using boost::ublas vectors is about 20% slower than Ublass non-generic version. However, if I substitute the container for a std::array, the generic-algorithm is 40% faster.
Below, listing 1. is an implementation of a generic algorithm – the Runga-kutta method which is well known for its use in the approximation of solutions of ordinary differential equations. Although the Runga-kutta method is well known and relatively simple, the generic algorithm implementation below is difficult to read. I won't attempt to explain how it works; I wrote it several months ago and I don't remember all the details. The reason why I list it here is to be able to compare it with a non-generic implementation of the algorithm I wrote in listing 2.
Listing 2. is far easier to read and understand. It is almost a carbon copy of the mathematical equations for the Runga-kutta4 method (see listing 3). The non-generic implementation is easy to read (assuming you are familiar with the equations) and it was very easy to write. I was done in about 30 minutes, whereas the generic version took me a couple of days to complete. So why make the effort ? The generic version allows you to use any container that supports forward iterators e.g. std::array, std::vector, C array, boost::ublas vectors, Blitz++. With the non-generic version you are tied into using only boost::ublas vectors.
As far as performance is concerned, the generic version using boost::ublas vectors is about 20% slower than Ublass non-generic version. However, if I substitute the container for a std::array, the generic-algorithm is 40% faster.
Runga-kutta4 Performance of the generic-algorithm using different containers compared to the non-generic version:
std::vector: 25% faster
Blitz: about the same
ublas: 20% slower
std::array: 40% faster
template <class InputIterator, class ForwardIterator, class T, class Time=double>
std::pair<std::vector<Time>, ForwardIterator> __runge_kutta4(
InputIterator inBegin, InputIterator inEnd,
Time x0, Time x1, int nr_steps, T f, ForwardIterator yout)
{
typedef typename std::iterator_traits<InputIterator>::value_type
value_type;
auto delta_t = static_cast<Time>((x1-x0)/nr_steps);
auto nr_equations = std::distance(inBegin, inEnd);
ForwardIterator youtP = yout;
yout = std::copy(inBegin, inEnd, yout);
std::vector<Time> ytime(nr_steps);
for (auto ii=0; ii< nr_steps; ++ii) {
typedef typename std::vector<value_type> vec;
// K1 = f(x, u(x))
vec v1(nr_equations);
typename vec::iterator k1 = f(yout-nr_equations, x0, begin(v1));
// K2 = f(x+deltaX/2, u(x)+K1*deltaX/2)
vec v2(nr_equations);
vec y2(nr_equations);
std::transform(yout-nr_equations, yout, k1, begin(y2),
[delta_t](value_type& y,value_type& k) {
return y+0.5*delta_t*k;
});
typename vec::iterator k2 =
f(begin(y2), 0.5*delta_t+x0, begin(v2));
// K3 = f(x+deltaX/2, u(x)+K2*deltaX/2)
vec v3(nr_equations);
vec y3(nr_equations);
std::transform(yout-nr_equations, yout, k2, begin(y3),
[delta_t](value_type& y,value_type& k){
return y+0.5*delta_t*k;
});
typename vec::iterator k3 =
f(begin(y3), 0.5*delta_t+x0, begin(v3));
// K4 = f(x+deltaX, u(x)+K2*deltaX)
vec v4(nr_equations);
vec y4(nr_equations);
std::transform(yout-nr_equations, yout, k3, begin(y4),
[delta_t](value_type& y, value_type& k){ return y+delta_t*k; });
typename vec::const_iterator k4
= f(begin(y4), delta_t+x0, begin(v4));
// u(x+deltaX) = u(x) + (1/6) (K1 + 2 K2 + 2 K3 + K4) * deltaX
for (auto i=0; i<nr_equations; ++i,++yout) {
*yout = *(yout-nr_equations) +
(1.0/6.0) * (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]) * delta_t;
}
ytime[ii] = x0;
x0 += delta_t;
}
return std::make_pair(ytime, yout);
}
Listing 1: Generic Implementation of Runga-kutta4
template <class T, class Vector=boost::numeric::ublas::vector<double> >
Vector rk4(const Vector& init, double x0, double x1, int nr_steps, T f)
{
Vector y = init;
auto h = static_cast<double>((x1-x0)/nr_steps);
double xn = x0;
for (auto ii=0; ii< nr_steps; ++ii) {
Vector k1 = f(y, xn);
Vector k2 = f(y+0.5*h*k1, 0.5*h+xn);
Vector k3 = f(y+0.5*h*k2, 0.5*h+xn);
Vector k4 = f(y+h*k3, h+xn);
xn += h;
y = u + (1.0/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4) * h;
}
return y;
}
Listing 2: Non-Generic Implementation of Runga-kutta4 (using ublas)
for n = 0,1,2,3,..., using
k1 = h f(xn, yn)
k2 = h f(xn + h/2, yn + k1/2)
k3 = h f(xn + h/2, yn + k2/2)
k4 = h f(xn + h, yn + k3)
xn+1 = xn + h
yn+1 = yn + h(1/6)(k1 + 2k2 + 2k3 + k4)
Listing 3: Runga-kutta4 equations.
วิดิโอสล็อตรีวิว POP เกมส์สล็อตทำเงิน ธีมของสวนสนุกที่เต็มไปด้วยลูกโป่ง ที่สามารถทำเงินได้อย่างหลายเท่าตัว โอกาสการทำเงินออนไลน์ ไม่ว่าเอเย่นแบบลื่นไหลผ่านทาง https://www.megaslot.game/ ทำเงินออนไลน์ กับสล็อค สิโนออนไลน์ 24 ชั่วโมง.
ReplyDelete