# Exact Pendulum modeling

1. Jan 28, 2004

### franznietzsche

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 $$x = x_0 + v*t + 1/2a*t^s$$ 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. Jan 28, 2004

### Integral

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

3. Jan 28, 2004

### nbo10

The pendulum problem can't be solved exactly. Aprroximation have to be used.

JMD

4. Jan 28, 2004

### krab

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 $x = x_0 + v*t$; neither is a good approximation. Look up Numerical Solutions of ODE's, especially Runge-Kutta methods.

Last edited: Jan 28, 2004
5. Jan 28, 2004

### franznietzsche

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:

$$d^2(theta)/(dt^2) = (g/l)*sin(theta)$$
which is the poalr coordinate form of the equations of motion.

6. Jan 28, 2004

### franznietzsche

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.

7. Jan 28, 2004

### Integral

Staff Emeritus
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.

8. Jan 28, 2004

### franznietzsche

the angular acceleration (listed as aac in the code) is :

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

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
9. Jan 28, 2004

### Integral

Staff Emeritus

10. Jan 29, 2004

### Integral

Staff Emeritus
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:

$$U(t)= \frac {d \theta} {dt}$$

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

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

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
11. Jan 29, 2004

### franznietzsche

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.

12. Jan 29, 2004

### Integral

Staff Emeritus
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
$$U_{n+1}= U_n - \frac {hg} {l} sin(\theta_n)$$

Note the change in sign for the last term. It should work now.

13. Jan 29, 2004

### nbo10

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
14. Jan 29, 2004

### nbo10

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

15. Jan 29, 2004

### franznietzsche

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

16. Jan 29, 2004

### franznietzsche

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?

17. Jan 29, 2004

### nbo10

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

18. Jan 29, 2004

### Integral

Staff Emeritus
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.

19. Jan 29, 2004

### franznietzsche

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.

20. Jan 29, 2004

### Integral

Staff Emeritus
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
$$U_{n+1}= U_n - \frac {hg} {l} sin(\theta_n)$$

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

Then
$$\theta_{n+1} = \theta_n + h U_{n+1}$$

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