Archived

This topic is now archived and is closed to further replies.

Runge Kutta - Oscillating Systems

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

Recommended Posts

Please, I want some help with the Runge Kutta method for solving differential equations of second order (y''''= -k y , k>0). I know there is an analytical solution for this, but I wanted to use Runge Kutta to simulate an oscillating system. The problem is that the Amplitude isn''t behaving constant as it should. I am using a small step, and the formula I am using is found in Hibbeler Physic Book. CAN SOMEONE PLEASE SEND ME THE METHOD IN C? I REALLY WANT YOUR HELP AND ANY ADVICE WOULD BE WELCOME. THANK YOU ALL IN ADVANCE, Daniel

Share on other sites
You''re correct, the oscillating system (let''s call it the spring equation from now on), is very sensitive to amplitude changes due to numerical error (which, even worse, cumulates after each step).

There are better numerical methods for the solution of the spring equation than Runge-Kutta 4 (I''m assuming you want the fourth order method). But anyway, here''s an untested source:

// this code is C++, as it makes this easy and clean
// the easiest way is to represent the "state" of the system
// is as a separate class - this way, we can use the "usual"
// first-order Runge-Kutta equations without modifications

typedef double real; // choose double or float
#define K 0.5 // choose your own value of K

class State
{
public:
State operator + ( const State & o )
{
State t;
t.y = y+o.y;
t.yDot = yDot+o.yDot;
return t;
}

State operator * ( const real a )
{
State t;
t.y = y*a;
t.yDot = yDot*a;
return t;
}

real y; // "position"
real yDot; // "velocity"
};

State Derivative( State x )
{
State t;
t.y = x.yDot; // derivative of position is velocity
t.yDot = -K*x.y; // derivative of velocity is acceleration
return t;
}

// returns the new state
// h = step size
State RK4( State x, real h )
{
State k1 = Derivative( x );
State k2 = Derivative( x+k1*(h/2.0) );
State k3 = Derivative( x+k2*(h/2.0) );
State k4 = Derivative( x+k3*h );

return x + (k1 + k2*2.0 + k3*2.0 + k4)*(h/6.0);
}

If you have any further questions, don''t hesitate to ask.

- Mikko Kauppila

1. 1
2. 2
Rutin
23
3. 3
4. 4
5. 5
khawk
14

• 9
• 11
• 11
• 23
• 12
• Forum Statistics

• Total Topics
633653
• Total Posts
3013160
×