/* 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);
}
*/
}
4th Order Runge Kutta
Author
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:
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.
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.
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.
Author
That's just it though. I want to use t, and it's not using it. What can I do with t? Should I multiply it with y?
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).
Author
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 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.
Author
Actually I planned on using this for Rigid Body physics, so it shouldn't be overkill. Never worked with rigid bodies before.
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.
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.
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.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
