Runge-Kutta Method for a double pendulum

[heycoa]

Hello, I am trying to program a double pendulum via the 4th order Runge-Kutta method and I cannot seem to be getting the right output. At first I used the Euler-Cromer method, but now I am aiming to make it more accurate.

Homework Statement

I have the equations of motion: d(omega)/dt and d(theta)/dt = omega

also, my step size is "h"

Homework Equations

The Runge-Kutta method can be found here: http://en.wikipedia.org/wiki/Runge–Kutta_methods

The Attempt at a Solution

for omega I tried the following:
k1=d(omega)/dt
k2=d(omega)/dt + 0.5*h*k1
k3=d(omega)/dt + 0.5*h*k2
k4=d(omega)/dt + 0.5*h*k3

omega= omega_initial + (1/6)*h*(k1 + 2*k2 + 2*k3 + k4)

and for theta:
k1=omega
k2=omega + 0.5*h*k1
k3=omega + 0.5*h*k2
k4=omega + 0.5*h*k3

theta = theta_initial + (1/6)*h*(k1 + 2*k2 + 2*k3 + k4)

 

[pasmith]

Hello, I am trying to program a double pendulum via the 4th order Runge-Kutta method and I cannot seem to be getting the right output. At first I used the Euler-Cromer method, but now I am aiming to make it more accurate.

Homework Statement

I have the equations of motion: d(omega)/dt and d(theta)/dt = omega

Homework Equations

The Runge-Kutta method can be found here: http://en.wikipedia.org/wiki/Runge–Kutta_methods

The Attempt at a Solution

for omega I tried the following:
k1=d(omega)/dt
k2=d(omega)/dt + 0.5*dt*k1
k3=d(omega)/dt + 0.5*dt*k2
k4=d(omega)/dt + 0.5*dt*k3

omega= omega_initial + (1/6)*dt*(k1 + 2*k2 + 2*k3 + k4)

and for theta:
k1=omega
k2=omega + 0.5*dt*k1
k3=omega + 0.5*dt*k2
k4=omega + 0.5*dt*k3

theta = theta_initial + (1/6)*dt*(k1 + 2*k2 + 2*k3 + k4)

Not quite; the k_i are vectors with one component for each variable, and you need to compute each component of k_1 before computing any components of k_2 and so forth. Thus if you have x = (x_1, \dots, x_N) and you wish to solve \dot x = f(x) then you would need (in C):

double k1[N], k2[N], k3[N], k4[N], xTemp[N], xNew[N];

for(i = 0; i < N; i++)
{
   k1[i] = f(xOld)[i];
}

for(i = 0; i < N; i++)
{
   xTemp[i] = xOld[i] + 0.5*dt*k1[i];
}

for(i = 0; i < N; i++)
{
   k2[i] = f(xTemp)[i]; 
}
for(i = 0; i < N; i++)
{
   xTemp[i] = xOld[i] + 0.5*dt*k2[i];
}

with similar loops to calculate k_3 and k_4.

Of course if you're using a language which understands vector operations then you don't need to expressly loop through the components.

 

[heycoa]

Thank you very much for the response! But I'm not sure I understand. What is N here?

 

[heycoa]

I guess I just don't understand what K is. In my case, is K the slope of d(omega)/dt? So do I have to take the derivative of d(omega)/dt and then calculate the slope at each point?