Jump to content
  • Advertisement
Sign in to follow this  
Omid Ghavami

Please criticize my RK4 implementation

This topic is 4449 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Hi! I just did an RK4 implementation, and I could use some help from you C++ gurus to improve the code somewhat :) I want the method interface to be as close to MATLAB's ode45 method as possible. Which means that I need variable dimensional vectors, and my own choice of solution to the problem is the one given below. Now I made the choice of using operators to make the method implementation code as clean as possible. However I know that it is not such a good practice (?). So I'm looking for advice on how to improve this code: Both in regard to good coding practice And also performance-wise. The solver method
MathVector Ode::Rk4(boost::function<MathVector (float t, MathVector& y)> f, MathVector& y, float t, float h)
{
	MathVector k1 = f(t, y);
	MathVector k2 = f(t + 0.5f * h, y + 0.5f * h * k1);
	MathVector k3 = f(t + 0.5f * h, y + 0.5f * h * k2);
	MathVector k4 = f(t + h, y + h * k3);

	return (y + (h / 6.0f) * (k1 + 2 * k2 + 2 * k3 + k4));
}

The MathVector header (std::vector override)
class MathVector :
	public std::vector<float>
{

	void operator += ( MathVector& v );
	void operator -= ( MathVector& v );

	friend MathVector operator + ( const MathVector& v1, const MathVector& v2 );
	friend MathVector operator - ( const MathVector& v1, const MathVector& v2 );
	friend MathVector operator * ( const MathVector& v, const float f );
	friend MathVector operator * ( const float f, const MathVector& v );
	friend MathVector operator * ( const MathVector& v1, const MathVector& v2 );
};

The MathVector implementation
/////////////////////////////// PUBLIC ///////////////////////////////////////
float add(const float a, const float b)
{
	return (a + b);
}

float sub(const float a, const float b)
{
	return (a - b);
}

float mul(const float a, const float b)
{
	return (a * b);
}

//============================= OPERATORS ====================================
void MathVector::operator += ( MathVector& v )
{
	std::transform(this->begin(), this->end(), v.begin(), this->begin(), boost::bind(&add, _1, _2));
}

void MathVector::operator -= ( MathVector& v )
{
	std::transform(this->begin(), this->end(), v.begin(), this->begin(), boost::bind(&sub, _1, _2));
}


MathVector operator + ( const MathVector& v1, const MathVector& v2 )
{
	MathVector result;
	result.reserve(v1.size());

	std::transform(v1.begin(), v1.end(), v2.begin(), std::back_inserter(result), boost::bind(&add, _1, _2));

	return result;
}

MathVector operator - ( const MathVector& v1, const MathVector& v2 )
{
	MathVector result;
	result.reserve(v1.size());

	std::transform(v1.begin(), v1.end(), v2.begin(), std::back_inserter(result), boost::bind(&sub, _1, _2));

	return result;
}

MathVector operator * ( const MathVector& v1, const float f )
{
	MathVector result;
	result.reserve(v1.size());

	std::transform(v1.begin(), v1.end(), std::back_inserter(result), boost::bind(&mul, _1, f));

	return result;
}

MathVector operator * ( const float f, const MathVector& v1 )
{
	MathVector result;
	result.reserve(v1.size());

	std::transform(v1.begin(), v1.end(), std::back_inserter(result), boost::bind(&mul, _1, f));

	return result;
}

MathVector operator * ( const MathVector& v1, const MathVector& v2 )
{
	MathVector result;
	result.reserve(v1.size());

	std::transform(v1.begin(), v1.end(), v2.begin(), std::back_inserter(result), boost::bind(&mul, _1, _2));

	return result;
}

If I missed posting any significant code, then just ask :) Any help appreciated! Thanks in advance

Share this post


Link to post
Share on other sites
Advertisement
I might just be delusional at this hour of the night, but you might want to look into boost::lambda instead of all those boost::bind things and functions like "add" :) though perhaps there's a reason for bind I'm not seeing?

Share this post


Link to post
Share on other sites
Your Ode::Rk4 function looks nice. But you might want to make it a template function so you can use it with any datatype that supports basic arithmetics. Your input vector (MathVector& y) should be declared const if you don't want to modify it.

Should MathVector be (publicly) derived from std::vector, or should it contain a std::vector as member ? Also you may want to make MathVector a template, having the storage type as template parameter. This way you can operate with vectors storing elements of any type that supports the necessary arithmetic operations.

Share this post


Link to post
Share on other sites
Hi again!

Sorry for the late reply, I got caught up in study :P

Thanks to both of you for the great advice! Rating++

etothex: Boost::Lambda seems like a much better solution indeed, I must admit it was my lack of boost knowledge which led to the usage of bind.

Quote:

But you might want to make it a template function so you can use it with any datatype that supports basic arithmetics


nmi: Templating the rk4 function sounds like a great idea. I am however unsure as to how to put a requirement such as "any datatype that supports basic arithmetics". Any help on that appreciated!

Also another question which I was wondering about. In my code I always left multiply scalar with object. Are there any reasons why either left or right multiplication would be prefered over the other?

Thanks in advance

Share this post


Link to post
Share on other sites
Quote:
Original post by Omid Ghavami
nmi: Templating the rk4 function sounds like a great idea. I am however unsure as to how to put a requirement such as "any datatype that supports basic arithmetics". Any help on that appreciated!



You can spedify template parameters for your function, which the compiler will try to deduce. For instance take this version of RK4:

template <typename F, typename T>
T Rk4(F& f, T& y, float t, float h)
{
T k1 = f(t, y);
T k2 = f(t + 0.5f * h, y + 0.5f * h * k1);
T k3 = f(t + 0.5f * h, y + 0.5f * h * k2);
T k4 = f(t + h, y + h * k3);

return (y + (h / 6.0f) * (k1 + 2 * k2 + 2 * k3 + k4));
}



Here two template parameters (F and T) were provided. The function expects a reference to something that is of type F as first parameter, and it expects as second parameter something of type T. (Maybe you also want to make the type of t and h a template parameter or deduce their type from a typedef provided by T.)

When you call the function, you have to provide values that allow the operations like addition (for y) or function call (for f):

Functor func;
MathVector v;
float time;
float step;

MathVector result = Rk4(func, v, time, step);



Instead of functor (template parameter F) you can use anything that allows an object of that type to be called like a function (i.e. it is a function or it is a class with overloaded "operator()" ). Same goes for template parameter T.

This way your Rk4 only describes the algorithm. It does not depend on any types.


Share this post


Link to post
Share on other sites
Quote:

class MathVector : public std::vector<float>


Oh dear. This is not a good idea.

First of all, std::vector does not provide a virtual destructor so if you ever happen to be in the (admittedly unlikely) situation of having to delete a MathVector* through a std::vector*, bad things will happen.

Second, std::vector - dispite the name -- really has nothing to do with a mathmetical vector. std::vector models a resizable array. By doing this you are allowing completely nonsensical behavior like "MathVector v; v.insert(whatever);"

A (mathematical) vector IS-GENERALLY-NOT a dynamic, resizeable array (though it could be implemented with one). Inheritance is not what you want here; you want your MathVector to contain a std::vector storing its elements (if even that; if your vectors are not actually meant to be re-dimensioned at runtime, a static array whose size is specified by a non-type template parameter might suffice).

It seems like you are going overboard on C++isms like boost::bind and such when you don't really need to. For example, all your add, sub, mul, et cetera functions have already been provided for you in a more robust fashion (std::plus, et al). So your implementations are redundant and sub-par (as they work only with floats).

Providing overloaded * for dot product (or ^ or anything else for cross product) is generally frowned upon (lack of clarity) although this is a stylistic choice.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!