1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Runge-Kutta Method for a double pendulum

  1. Nov 29, 2013 #1
    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.

    1. The problem statement, all variables and given/known data

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

    also, my step size is "h"

    2. Relevant equations

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

    3. 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)
     
    Last edited: Nov 29, 2013
  2. jcsd
  3. Nov 29, 2013 #2

    pasmith

    User Avatar
    Homework Helper

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

    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 [itex]k_3[/itex] and [itex]k_4[/itex].

    Of course if you're using a language which understands vector operations then you don't need to expressly loop through the components.
     
  4. Nov 29, 2013 #3
    Thank you very much for the response! But I'm not sure I understand. What is N here?
     
  5. Nov 29, 2013 #4
    I guess I just dont 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?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted