Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Exact Pendulum modeling

  1. Jan 28, 2004 #1
    I'm trying to use the differential equations for the motion of a pendulum to create a program that predcits the pendulums position.

    It uses the [tex] x = x_0 + v*t + 1/2a*t^s [/tex] and is written in C++.

    the code is:

    #include<iostream.h>
    #include<stdlib.h>
    #include<pendulumfunctions.hpp>

    int main()
    {
    cout << "Enter a value for g: ";
    cin >> g;
    cout << "\nEnter an initial angular position:";
    cin >> iang;
    ang = iang;
    y = ang;
    cout << "\nEnter the pendulum length: ";
    cin >> l;
    cout << "\n";
    system("PAUSE");
    t = 0;
    av = 0;
    while (t<=10)
    {
    ang = y;
    av = w;
    cout << "Time: " << t << "\tVelocity: " << velocity(g, l, ang) << "\tPosition: " << position(g, l, ang) << endl;
    t = t + (10.0/2500.0);
    }
    system("PAUSE");
    return 0;
    }


    PendulumFunctions.hpp :

    #include<math.h>

    typedef double dou;

    dou g, l, ang, iang, av, C, v, y, w, t;


    dou aacceleration(dou g, dou l, dou ang) //Function for angular acceleration, radians/sec
    {
    dou aa;
    aa = g/l*sin(ang);
    return (aa);
    }

    dou avelocity(dou g, dou l, dou ang)
    {
    if (t==0)
    return av;
    else
    {
    w = av + (10.0/2500.0)*aacceleration(g, l, ang);
    return (av);
    }
    }

    dou velocity(dou g, dou l, dou ang)
    {
    v = avelocity(g, l, ang)*l;
    return (v);
    }

    dou position(dou g, dou l, dou ang)
    {
    if (t==0)
    return iang;
    else
    {
    y = ang + 10.0/2500.0*avelocity(g, l, ang) + 0.5*(10.0/2500.0)*(10.0/2500.0)*aacceleration(g, l, ang);
    return (y);
    }
    }


    Where (10.0/2500.0) is the time increment from point to point, ang = angle, etc.

    NOw my question is when i run this, the answer starts out cycling as normal, but then the maximum ampltude and period sloy increase until at the final postion t = 10 seconds, the pendulum has swung all the way out to -18 radians, something that makes no sense for a conservative system driven only by gravity. Can anyone help with this?
     
  2. jcsd
  3. Jan 28, 2004 #2

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    The equation of motion you have used is not that of a pendulum. It is equation for a free falling body, if a=g.

    What is your Differential Equation?
     
  4. Jan 28, 2004 #3
    The pendulum problem can't be solved exactly. Aprroximation have to be used.

    JMD
     
  5. Jan 28, 2004 #4

    krab

    User Avatar
    Science Advisor

    I don't know how you got any oscillatory behaviour, since you should have

    aa= -g/l*sin(ang);

    (Acceleration is towards ang=0! Edit: Alternatively, you input g as a negative number.) As you stated, you then get oscillatory behaviour of slowly growing amplitude. I get 10% growth after 5 periods. The reason is that you are using a constant acceleration in each time step. You could as well use a constant velocity; simply [itex] x = x_0 + v*t[/itex]; neither is a good approximation. Look up Numerical Solutions of ODE's, especially Runge-Kutta methods.
     
    Last edited: Jan 28, 2004
  6. Jan 28, 2004 #5
    First g is entered as a negative, so the negative sign is built into g, not the equation. Also, the acceleration is only constant in each iteration. the value of ang changes from iteration to iteration (becoming the y value). Thanx for pointing me to the Runge-Kutta method, hven't heard of it before, but i will check it out.


    I know it can't belsoved exactly, this is the method of approximation that i am trying to figure out. Its not an exact solution.

    I'm using polar coordinates for the differential equation in which case it is a=(g/l)sin(theta). Even if i wasn't a=g is wrong because it does not take into acount the restrained motion caused by the force of the tension in the cord. the equation of motion and differential equation are one and the same because:

    [tex] d^2(theta)/(dt^2) = (g/l)*sin(theta) [/tex]
    which is the poalr coordinate form of the equations of motion.
     
  7. Jan 28, 2004 #6
    I took your advice and used the Runge-Kutta Algorithm. However, for some reason my results went from having oscillatory behavoir to just being outright divergent:

    #include<iostream.h>
    #include<stdlib.h>
    #include<pendulumfunctions.h>

    int main()
    {
    cout << "Enter a value for g: ";
    cin >> g;
    cout << "\nEnter an initial angular position:";
    cin >> ang;
    iang = ang;
    y = ang;
    cout << "\nEnter the pendulum length: ";
    cin >> l;
    cout << "\n";
    system("PAUSE");
    t = 0;
    av = 0;
    w = 0;
    h = .1;
    while (t<=10)
    {
    ang = y;
    av = w;
    cout << "Time: " << t << "\tVelocity: " << velocity(g, l, ang) << "\tPosition: " << position(g, l, ang) << endl;
    t = t + (h);
    }
    system("PAUSE");
    return 0;
    }


    PENDULUMFUNCTIONS.h:

    #include<math.h>

    typedef double dou;

    dou g, l, ang, iang, av, aa, t, h, w, y, v;


    dou aac(dou g, dou l, dou ang) //Function for angular acceleration, radians/sec
    {
    dou aa;
    aa = g/l*sin(ang);
    return (aa);
    }


    dou vkone()
    {
    return(aac(g, l, ang));
    }

    dou vktwo()
    {
    return(aac(g, l, (ang + (h/2)*vkone())));
    }

    dou vkthree()
    {
    return(aac(g, l, (ang + (h/2)*vktwo())));
    }

    dou vkfour()
    {
    return(aac(g, l, (ang + h*vkthree())));
    }

    dou ave(dou g, dou l, dou ang) //Uses Runge-Kutta Algorithm for angular velocity
    {
    if (t==0)
    return av;
    else
    {
    w = ang +(h/6)*(vkone() + 2*vktwo() + 2*vkthree() + vkfour());
    return w;
    }
    }

    dou velocity(dou g, dou l, dou ang)
    {
    v = ave(g, l, ang)*l;
    return (v);
    }

    dou pkone()
    {
    return(ave(g, l, ang));
    }

    dou pktwo()
    {
    return(ave(g, l, ang + (h/2)*pkone()));
    }

    dou pkthree()
    {
    return(ave(g, l, ang + (h/2)*pktwo()));
    }

    dou pkfour()
    {
    return(ave(g, l, (ang + h*pkthree())));
    }

    dou position(dou g, dou l, dou ang)
    {
    if (t==0)
    return iang;
    else
    {
    y = ang + (h/6)*(pkone() + 2*pktwo() + 2*pkthree() + pkfour());
    return y;
    }
    }


    OUTPUT:

    Enter a value for g: -9.8

    Enter an initial angular position:0.785398

    Enter the pendulum length: 1

    Press any key to continue . . .
    Time: 0 Velocity: 0 Position: 0.785398
    Time: 0.1 Velocity: 0.312324 Position: 0.818245
    Time: 0.2 Velocity: 0.326553 Position: 0.852589
    Time: 0.3 Velocity: 0.341597 Position: 0.888515
    Time: 0.4 Velocity: 0.357526 Position: 0.926117
    Time: 0.5 Velocity: 0.374424 Position: 0.965495
    Time: 0.6 Velocity: 0.392384 Position: 1.00676
    Time: 0.7 Velocity: 0.411514 Position: 1.05004
    Time: 0.8 Velocity: 0.431942 Position: 1.09547
    Time: 0.9 Velocity: 0.453816 Position: 1.1432
    Time: 1 Velocity: 0.477313 Position: 1.1934
    Time: 1.1 Velocity: 0.502645 Position: 1.24626
    Time: 1.2 Velocity: 0.530065 Position: 1.30201
    Time: 1.3 Velocity: 0.559885 Position: 1.36089
    Time: 1.4 Velocity: 0.592489 Position: 1.4232
    Time: 1.5 Velocity: 0.628356 Position: 1.48929
    Time: 1.6 Velocity: 0.668092 Position: 1.55955
    Time: 1.7 Velocity: 0.712476 Position: 1.63448
    Time: 1.8 Velocity: 0.762522 Position: 1.71468
    Time: 1.9 Velocity: 0.819576 Position: 1.80087
    Time: 2 Velocity: 0.885461 Position: 1.894
    Time: 2.1 Velocity: 0.962704 Position: 1.99525
    Time: 2.2 Velocity: 1.0549 Position: 2.10619
    Time: 2.3 Velocity: 1.16735 Position: 2.22896
    Time: 2.4 Velocity: 1.30813 Position: 2.36654
    Time: 2.5 Velocity: 1.49013 Position: 2.52326
    Time: 2.6 Velocity: 1.73479 Position: 2.70571
    Time: 2.7 Velocity: 2.07845 Position: 2.9243
    Time: 2.8 Velocity: 2.57773 Position: 3.1954
    Time: 2.9 Velocity: 3.28428 Position: 3.54081
    Time: 3 Velocity: 4.12749 Position: 3.9749
    Time: 3.1 Velocity: 4.87354 Position: 4.48746
    Time: 3.2 Velocity: 5.4018 Position: 5.05557
    Time: 3.3 Velocity: 5.76288 Position: 5.66166
    Time: 3.4 Velocity: 6.03982 Position: 6.29687
    Time: 3.5 Velocity: 6.28841 Position: 6.95823
    Time: 3.6 Velocity: 6.54873 Position: 7.64696
    Time: 3.7 Velocity: 6.8773 Position: 8.37026
    Time: 3.8 Velocity: 7.43022 Position: 9.1517
    Time: 3.9 Velocity: 8.7251 Position: 10.0693
    Time: 4 Velocity: 10.8756 Position: 11.2131
    Time: 4.1 Velocity: 11.9782 Position: 12.4729
    Time: 4.2 Velocity: 12.5307 Position: 13.7907
    Time: 4.3 Velocity: 13.085 Position: 15.1669
    Time: 4.4 Velocity: 14.4385 Position: 16.6854
    Time: 4.5 Velocity: 17.6181 Position: 18.5383
    Time: 4.6 Velocity: 18.73 Position: 20.5082
    Time: 4.7 Velocity: 19.6288 Position: 22.5725
    Time: 4.8 Velocity: 23.3337 Position: 25.0266
    Time: 4.9 Velocity: 25.0922 Position: 27.6655
    Time: 5 Velocity: 26.8838 Position: 30.4929
    Time: 5.1 Velocity: 31.0429 Position: 33.7577
    Time: 5.2 Velocity: 32.8712 Position: 37.2148
    Time: 5.3 Velocity: 37.5114 Position: 41.1599
    Time: 5.4 Velocity: 41.6486 Position: 45.5401
    Time: 5.5 Velocity: 44.6937 Position: 50.2406
    Time: 5.6 Velocity: 50.256 Position: 55.5261
    Time: 5.7 Velocity: 56.1297 Position: 61.4293
    Time: 5.8 Velocity: 62.2155 Position: 67.9726
    Time: 5.9 Velocity: 68.6381 Position: 75.1913
    Time: 6 Velocity: 75.319 Position: 83.1126
    Time: 6.1 Velocity: 82.3145 Position: 91.7697
    Time: 6.2 Velocity: 92.588 Position: 101.507
    Time: 6.3 Velocity: 100.928 Position: 112.122
    Time: 6.4 Velocity: 112.7 Position: 123.975
    Time: 6.5 Velocity: 124.863 Position: 137.107
    Time: 6.6 Velocity: 137.763 Position: 151.595
    Time: 6.7 Velocity: 151.115 Position: 167.488
    Time: 6.8 Velocity: 168.422 Position: 185.201
    Time: 6.9 Velocity: 184.953 Position: 204.653
    Time: 7 Velocity: 205.294 Position: 226.244
    Time: 7.1 Velocity: 226.213 Position: 250.035
    Time: 7.2 Velocity: 250.773 Position: 276.409
    Time: 7.3 Velocity: 276.441 Position: 305.482
    Time: 7.4 Velocity: 306.347 Position: 337.701
    Time: 7.5 Velocity: 338.559 Position: 373.308
    Time: 7.6 Velocity: 372.579 Position: 412.492
    Time: 7.7 Velocity: 413.419 Position: 455.972
    Time: 7.8 Velocity: 456.604 Position: 503.993
    Time: 7.9 Velocity: 503.235 Position: 556.919
    Time: 8 Velocity: 557.825 Position: 615.586
    Time: 8.1 Velocity: 615.688 Position: 680.338
    Time: 8.2 Velocity: 679.433 Position: 751.795
    Time: 8.3 Velocity: 752.724 Position: 830.959
    Time: 8.4 Velocity: 830.106 Position: 918.262
    Time: 8.5 Velocity: 917.715 Position: 1014.78
    Time: 8.6 Velocity: 1014.85 Position: 1121.51
    Time: 8.7 Velocity: 1121.45 Position: 1239.46
    Time: 8.8 Velocity: 1238.57 Position: 1369.72
    Time: 8.9 Velocity: 1369.73 Position: 1513.77
    Time: 9 Velocity: 1514.06 Position: 1673.01
    Time: 9.1 Velocity: 1672.12 Position: 1848.87
    Time: 9.2 Velocity: 1848 Position: 2043.22
    Time: 9.3 Velocity: 2042.54 Position: 2258.04
    Time: 9.4 Velocity: 2257.17 Position: 2495.43
    Time: 9.5 Velocity: 2494.83 Position: 2757.81
    Time: 9.6 Velocity: 2758.12 Position: 3047.88
    Time: 9.7 Velocity: 3047.55 Position: 3368.4
    Time: 9.8 Velocity: 3368.03 Position: 3722.62
    Time: 9.9 Velocity: 3722.34 Position: 4114.1
    Time: 10 Velocity: 4114.88 Position: 4546.86
    Press any key to continue . . .

    These results clearly are incorrect, however as best i can tell it seems that i followed the alogrithm correctly. Does anyone know what i programmed incorrectly? Is it my implementation of the algorithm or something else? Thanx in advance.
     
  8. Jan 28, 2004 #7

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    I find it rather tedious to sort through your code. I must be looking for coding errors while I am attempting to sort out exactly what your equations of motion are. Could you please post first the basic equations you are approximating then outline (flow chart?) your computations. Once we work out that your starting place is correct we can start finding the error.
     
  9. Jan 28, 2004 #8
    the angular acceleration (listed as aac in the code) is :

    [tex]\frac {d^2 \theta} {dt^2} = (g/l)*sin(\theta) [/tex]

    the others follow from this nonlinear second order differential equation.

    The pkone-pkfour and vkone-vkfour are the k values in the Runge-Kutta algorithm, vk0ne-four being used to calculate the angular velocity from the angular acceleration, pkone-pkfour to calculate the position from that angular velocity value. (its all in the pendulumfunctions.h section).
     
    Last edited by a moderator: Jan 28, 2004
  10. Jan 28, 2004 #9

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    Please tell me what "follows".
     
  11. Jan 29, 2004 #10

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    Since you cannot tell me what you are doing, I must assume that you do not know what you are doing. Drop back to the Euler method, you computations failed,not because of the method but because you simply were not modeling correctly.

    Both Eulers and Runga Cutta are FIRST order methods. You cannot directly solve a 2nd order equation using them. The first step is to reduce your 2nd order equation to a system of First order equations. In this case if you let:

    [tex] U(t)= \frac {d \theta} {dt} [/tex]

    You can express the original DQ as
    [tex] \frac {d U(t)} {dt} = - \frac g l sin(\theta) [/tex]

    Now using a simple Eulers method you arrive at the following eqations to be computed
    [tex] \theta_{n+1} = \theta_n + h U_n [/tex]
    [tex] U_{n+1}= U_n - \frac {hg} {l} sin(\theta_n)[/tex]


    here is a table of values I computed for h= .1 l=1f, with intial conditions Theata(0)=.1 U(0)=0

    The first column is time, 2nd velocity, 3rd position

    0 0 0.1
    0.1 -0.097836748 0.090216325
    0.2 -0.186128865 0.071603439
    0.3 -0.256240288 0.04597941
    0.4 -0.301284235 0.015850986
    0.5 -0.316817551 -0.015830769
    0.6 -0.301304046 -0.045961173
    0.7 -0.256277952 -0.071588968
    0.8 -0.186180673 -0.090207036
    0.9 -0.097897623 -0.099996798
    1 -6.39968E-05 -0.100003198
    1.1 0.09777587 -0.090225611
    1.2 0.186077049 -0.071617906
    1.3 0.256202614 -0.045997644
    1.4 0.301264412 -0.015871203
    1.5 0.316817538 0.01581055
    1.6 0.301323844 0.045942935
    1.7 0.256315605 0.071574495
    1.8 0.186232474 0.090197743
    1.9 0.097958494 0.099993592
    2 0.000127994 0.100006392
    2.1 -0.097714987 0.090234893



    It is always a good idea to start with a set of easily verifialble numbers. I can see that inspite of using a simple Euler method and a pretty agressive step size the method is yielding belivable numbers.

    Edit:
    Fixed a sign error.
     
    Last edited: Jan 29, 2004
  12. Jan 29, 2004 #11
    thanx, that should fix the problem.

    Did you compute this by hand or by program? if so could you post the code? I tried running it myself, but it diverged rapidly giving -46 radians at 10 seconds.
     
  13. Jan 29, 2004 #12

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    This is a cut and paste from Excel. Unfortunatly, I did this at work last night and do not have the spread sheet or my notes in front of me. But I think I see a typo in my last post, the Velocity equation should be
    [tex] U_{n+1}= U_n - \frac {hg} {l} sin(\theta_n)[/tex]


    Note the change in sign for the last term. It should work now.
     
  14. Jan 29, 2004 #13
    You can type code in a code mode with [{code}] .. [{ / code}]. Take out the curly brakets.
    Code (Text):

    Put code here and it's much easier to read. The spacing and indents are much better.

     
    for future reference.

    JMD
     
    Last edited: Jan 29, 2004
  15. Jan 29, 2004 #14
    Code (Text):

    calculate(theta(),omega(),t(),length,dt){
       i = 0
       g = 9.8
       period = 2 * pi / sqr(g/length)    ! period of pendulum
       do
          i = i + 1
          t(i+1) = t(i) + dt
          omega(i+1) = omega(i) - (g/length) * theta(i)
          theta(i+1) = theta(i) + omega(i) * dt
       loop until t(i+1) >= 5 * period
    }
     
    JMD
     
  16. Jan 29, 2004 #15
    Using h=0.02 i get from a staring position of 0.785398 radians i get a position of 1.71817 radians at 8.88 seconds, which is over 2*PI radians
     
  17. Jan 29, 2004 #16
    Is that FORTRAN? i don't recognize the language.

    I tried to comprehend it, but i'm assuming omega is the angular velocity, but that doesn't seem to make sense. COuld you please explain?
     
  18. Jan 29, 2004 #17
    Opps I let out a dt
    Code (Text):

    calculate(theta(),omega(),t(),length,dt){
       i = 0
       g = 9.8
       period = 2 * pi / sqr(g/length)    ! period of pendulum
       do
          i = i + 1
          t(i+1) = t(i) + dt
          omega(i+1) = omega(i) - (g/length) * theta(i) * dt
          theta(i+1) = theta(i) + omega(i) * dt
       loop until t(i+1) >= 5 * period
    }

     
    it's basic, but I use it for more of a puesdocode.
    Now it should be clear

    JMD
     
  19. Jan 29, 2004 #18

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    There is no sense in computing the period of the pendulum. The equations I have provided will give a reasonable approximation to the pendulums motion for ALL beginning positions. The period you are computing results from a small angle approximation and is only valid for starting positions less then about .25 Radians. You should be able to see the difference in your computations.


    Frans,
    From what I can puzzle out of your code you are simply not using the Rk method correctly. You MUST start with the 2 DEQs which I posted above. Use the numbers generated to step forward to your next position. Once again I will ask you to show us the equations which you are being approximated by the RK (or Euler) methods. The more information you provide the easier it will be to help you.
     
  20. Jan 29, 2004 #19
    Well i'mt trying to create a computer model for a Group 4 lab project in my IB physics class. So no there is no real practical purpose or sense, but its still a problem to be soved.

    For my equations i use angular acceleration = (g/l)*sin(theta), then i use the Runge Kutta Algorithm on that to find the value of v the angular velocity, and then i use the Runge Kutta algorithm on those values to compute the angular position. Thats all in the section under pendulumfunctions.h. So i use the RK algorith to calculate velocity from acceleration, and position from my velocity result.

    Thanx for telling me about the restiction on the angle, i wil try using that and see how it turns out. Also, i tried your euler's method but it gave the same problem, diverging to about 11 radians after ten seconds.
     
  21. Jan 29, 2004 #20

    Integral

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    If you are still using the posted code for your RK approximation you will indeed get garbage. It simply is not correct. Drop the RK, go back to the Euler method. Once you figure out how to make it work you can then attempt the more involved RK computation.

    I use
    [tex] U_{n+1}= U_n - \frac {hg} {l} sin(\theta_n)[/tex]


    To generate a new velocity which is based upon the old position and velocity.

    Then
    [tex] \theta_{n+1} = \theta_n + h U_{n+1} [/tex]


    To generate a new postion based on the new velocity and old position.

    You cannot successfully use any numerical method by closing your eyes and cranking a handle. You must understand every computation you do and why you do it. I do not feel that you are at that point.

    EDIT: I have made major changes to the above post based up double checking my spread sheet. The above is the computations used to generate the numbers I posted last night.
     
    Last edited: Jan 29, 2004
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Exact Pendulum modeling
  1. Pendulum ? (Replies: 4)

  2. Pendulum's equilibrium (Replies: 3)

Loading...