Jump to content
  • Advertisement
Sign in to follow this  
Jacob Roman

4th Order Runge Kutta

This topic is 4651 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

I found some simple C code for Runge Kutta from here: http://www.physics.orst.edu/~rubin/nacphy/ComPhys/DIFFEQ/PRO/srk_c/srk_c.c But I feel it's wrong. Mainly cause f(x, y) isn't doing anything with the time. All it returns is a negative y. I converted most of it by change x into t for time, and got rid of saving it to file, but need some proper code for the f(t,y), or any other corrections that may be needed:
/* Runge Kutta algorithm for first-order differential equations*/

#include <stdio.h>

#define dist 0.5                        /* stepsize in t */ 
#define MAX 4                           /* max for t */

double runge4(double t, double y, double step); /* Runge-Kutta function */

double f(double t, double y);           /* function for derivative */

double  f(double t, double y)
{
	return(-y);
}

double runge4(double t, double y, double step)
{
	double h=step/2.0, k1, k2, k3, k4;

	k1=step*f(t, y);
	k2=step*f(t+h, y+k1/2.0);
	k3=step*f(t+h, y+k2/2.0);
	k4=step*f(t+step, y+k3);

	printf("%f\n",k1);
	printf("%f\n",k2);
	printf("%f\n",k3);
	printf("%f\n",k4);

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


main()
{
	double t, y;
	int j;

	y=1;                                    /* initial position */

	//printf("0\t%f\n", y);

	y = runge4(1, y, 10);

	/*
		for (j=1; j*dist<=MAX; j++)
		{
		t=j*dist;
		y=runge4(t, y, dist);

		printf("%f\t%f\n", t, y);
		}
	*/
}

Share this post


Link to post
Share on other sites
Advertisement
Guest Anonymous Poster
It isn't wrong. It's using the runge-kutta formula to integrate dy/dx = -y.

The output should be an exponential curve (it should exponentially decrease towards 0). If it isn't, then something is wrong.

You need to replace f with whatever you're integrating.

Share this post


Link to post
Share on other sites
y' = f(t, y) is the general form for an ODE, but many common applications--especially dynamics--won't use t at all. Consider, for instance, an equation describing the acceleration of a pendulum. The acceleration will vary given the bob's position and velocity, but not based on whether it's 3:00 or 5:00.

Share this post


Link to post
Share on other sites
Well, I don't know. Should your f(t,y) be t+y? Should it be 4t^2+7y? This has nothing to do with RK4 and everything to do with what it is you are actually simulating. f(t,y) is a function YOU will write, which represents your system. f(t,y) should be the derivative dy/dt at the point (t,y).

Share this post


Link to post
Share on other sites
Ok, well let's say I wanted y = y0 + a * t, which is a basic physics equation. Would it be ok for the function f(t,y) to return y0 + a * t, assuming I had y0 as an arguement?

Share this post


Link to post
Share on other sites
Quote:
Original post by Jacob Roman
Ok, well let's say I wanted y = y0 + a * t, which is a basic physics equation. Would it be ok for the function f(t,y) to return y0 + a * t, assuming I had y0 as an arguement?


Quote:
Original post by Sneftel
f(t,y) should be the derivative dy/dt at the point (t,y).


ergo, your function f should simply return a.

and yes, since 'a' is a constant, 4-th order runge kutta is overkill for this simple problem.

Share this post


Link to post
Share on other sites
Just to make sure that the issue is clear...

A set of first-order ordinary differential equations (ODE) can be written in explicit form as such:

(1) dy/dt = f(t, y),
where y is a vector and f is a vector function.

This means that dy/dt (or y') can be given as a function of t and y.

For example, let's look at Newton's equation of motion.
If the force acting on a point mass can be expressed as a function of the points position, we can write the equation of motion like this:

(2) dx/dt = v
(3) dv/dt = m-1F(x).

Now, if we define:
y = (x, v) and f(t, y) = (v, m-1F(x)),

then (1) is equivalent to (2) and (3).

For instance, if we have a satellite orbiting a planet which is located at the origin, we can write F(x) like this:
F(x) = -GMmx / |x|3

Here's a more empirical example:
(4) y(t) = t*e2t

The derivative of y w.r.t t is:
(5) dy/dt = e2t + 2*t*e2t = 1/t*y(t) + 2*y(t)

So if we define:
f(t,y) = 1/t*y(t) + 2*y(t),

we can write (5) like this:
(6) dy/dt = f(t,y).

One of the solutions to the ODE (6) is the equation (4).

PS. Most commerical engines use semi-implicit integration for rigid body simulation.

Share this post


Link to post
Share on other sites
Quote:
Original post by Jacob Roman
Actually I planned on using this for Rigid Body physics, so it shouldn't be overkill. Never worked with rigid bodies before.


rigid body physics can mean a lot of things, but regardless, i doubt runge kutta is the best solution for any rigid body problem i can imagine.

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!