Predicting Pendulum Position: Differential Equations in C++ Code

[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?

 

[Integral]

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?

 

[nbo10]

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

JMD

 

[krab]

Originally posted by 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^2$$ and is written in C++.

the code is:

#include<iostream.h> ...

dou aa;
aa = g/l*sin(ang); ...

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.

 

[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.

 

[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.

 

[Integral]

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.

 

[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).

 

[Integral]

Please tell me what "follows".

 

[Integral]

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 aggressive step size the method is yielding belivable numbers.

Edit:
Fixed a sign error.

 

[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.

 

[Integral]

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.

 

[nbo10]

You can type code in a code mode with [{code}] .. [{ / code}]. Take out the curly brakets.

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

for future reference.

JMD

 

[nbo10]

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

 

[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

 

[franznietzsche]

Originally posted by nbo10

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

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?

 

[nbo10]

Opps I let out a dt

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

 

[Integral]

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.

 

[franznietzsche]

Originally posted by Integral
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.

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.

 

[Integral]

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.

 

[Integral]

A side note concerning this type of calculation:
i = i + 1
t(i+1) = t(i) + dt

This is a recipe for disaster, your t(n) will accumulate the round off error inherent in dt with every pass. On my old 8 bit AppleII it only took about 10 passes for this error to show up. Recall that .1 cannot be stored precisely in a digital computer, so if you like to take steps in powers of .1 then you are creating problems for yourself. It would be better to take steps which are a power of .5 as they have no round off error.

If you wish to keep a running sum as above it is better to compute

i= i+1
t(i)= i * dt

This will give the same results while retaining accuracy.

And yes these are still meaningful concerns even with modern high soot computers.

 

[Integral]

There is an error in your RK code here

w = ang +(h/6)*(vkone() + 2*vktwo() + 2*vkthree() + vkfour());

This is your velocity computation, but you add the correction to the POSITION! you need to replace ang with the last computed velocity.

This error could cause some significant troubles.

 

[franznietzsche]

Thank you very much for ctaching that, it actually solved the entire problem, and the program gives very consistent data now. Thank you for all the help.

 

[krab]

Aww gee, why don't you post the output like you did before? Would be interesting..

 

[Integral]

I would recommend that you set g as a constant rather then as an input value, then compute g/l once, then pass that value rather then doing the same computation hundreds of times. Also you might consider making your initial velocity an input rather then setting it to 0 in the code.

 

[franznietzsche]

Originally posted by krab
Aww gee, why don't you post the output like you did before? Would be interesting..

Output:

Enter a value for g: -9.8

Enter an initial angular position:0.785398

Enter the pendulum length: 1

Enter the time interval between predictions: 0.05

Enter run-time length of simulation: 5

Press any key to continue . . .
Time: 0 Velocity: 0 Position: 0.785398
Time: 0.05 Velocity: -0.288253 Position: 0.770985
Time: 0.1 Velocity: -0.571795 Position: 0.742396
Time: 0.15 Velocity: -0.845896 Position: 0.700101
Time: 0.2 Velocity: -1.1058 Position: 0.644811
Time: 0.25 Velocity: -1.34676 Position: 0.577473
Time: 0.3 Velocity: -1.56412 Position: 0.499266
Time: 0.35 Velocity: -1.75342 Position: 0.411595
Time: 0.4 Velocity: -1.91054 Position: 0.316069
Time: 0.45 Velocity: -2.0319 Position: 0.214474
Time: 0.5 Velocity: -2.11462 Position: 0.108743
Time: 0.55 Velocity: -2.15668 Position: 0.000908989
Time: 0.6 Velocity: -2.15703 Position: -0.106942
Time: 0.65 Velocity: -2.11566 Position: -0.212726
Time: 0.7 Velocity: -2.03361 Position: -0.314406
Time: 0.75 Velocity: -1.91288 Position: -0.41005
Time: 0.8 Velocity: -1.75634 Position: -0.497867
Time: 0.85 Velocity: -1.56755 Position: -0.576245
Time: 0.9 Velocity: -1.35062 Position: -0.643776
Time: 0.95 Velocity: -1.11002 Position: -0.699277
Time: 1 Velocity: -0.850394 Position: -0.741796
Time: 1.05 Velocity: -0.576492 Position: -0.770621
Time: 1.1 Velocity: -0.293069 Position: -0.785275
Time: 1.15 Velocity: -0.00485676 Position: -0.785517
Time: 1.2 Velocity: 0.283435 Position: -0.771346
Time: 1.25 Velocity: 0.567096 Position: -0.742991
Time: 1.3 Velocity: 0.841394 Position: -0.700921
Time: 1.35 Velocity: 1.10158 Position: -0.645842
Time: 1.4 Velocity: 1.3429 Position: -0.578697
Time: 1.45 Velocity: 1.56069 Position: -0.500663
Time: 1.5 Velocity: 1.75049 Position: -0.413138
Time: 1.55 Velocity: 1.90818 Position: -0.317729
Time: 1.6 Velocity: 2.03017 Position: -0.216221
Time: 1.65 Velocity: 2.11356 Position: -0.110543
Time: 1.7 Velocity: 2.15631 Position: -0.00272696
Time: 1.75 Velocity: 2.15737 Position: 0.105141
Time: 1.8 Velocity: 2.1167 Position: 0.210976
Time: 1.85 Velocity: 2.03532 Position: 0.312742
Time: 1.9 Velocity: 1.91521 Position: 0.408503
Time: 1.95 Velocity: 1.75924 Position: 0.496465
Time: 2 Velocity: 1.57097 Position: 0.575014
Time: 2.05 Velocity: 1.35448 Position: 0.642737
Time: 2.1 Velocity: 1.11423 Position: 0.698449
Time: 2.15 Velocity: 0.854887 Position: 0.741193
Time: 2.2 Velocity: 0.581187 Position: 0.770253
Time: 2.25 Velocity: 0.297885 Position: 0.785147
Time: 2.3 Velocity: 0.0097135 Position: 0.785633
Time: 2.35 Velocity: -0.278616 Position: 0.771702
Time: 2.4 Velocity: -0.562393 Position: 0.743582
Time: 2.45 Velocity: -0.836888 Position: 0.701738
Time: 2.5 Velocity: -1.09735 Position: 0.64687
Time: 2.55 Velocity: -1.33902 Position: 0.579919
Time: 2.6 Velocity: -1.55725 Position: 0.502057
Time: 2.65 Velocity: -1.74756 Position: 0.414679
Time: 2.7 Velocity: -1.90582 Position: 0.319388
Time: 2.75 Velocity: -2.02843 Position: 0.217966
Time: 2.8 Velocity: -2.11249 Position: 0.112342
Time: 2.85 Velocity: -2.15594 Position: 0.00454492
Time: 2.9 Velocity: -2.1577 Position: -0.10334
Time: 2.95 Velocity: -2.11772 Position: -0.209226
Time: 3 Velocity: -2.03701 Position: -0.311077
Time: 3.05 Velocity: -1.91754 Position: -0.406954
Time: 3.1 Velocity: -1.76214 Position: -0.495061
Time: 3.15 Velocity: -1.57437 Position: -0.573779
Time: 3.2 Velocity: -1.35832 Position: -0.641696
Time: 3.25 Velocity: -1.11844 Position: -0.697618
Time: 3.3 Velocity: -0.859377 Position: -0.740587
Time: 3.35 Velocity: -0.585878 Position: -0.76988
Time: 3.4 Velocity: -0.302698 Position: -0.785015
Time: 3.45 Velocity: -0.0145702 Position: -0.785744
Time: 3.5 Velocity: 0.273795 Position: -0.772054
Time: 3.55 Velocity: 0.557688 Position: -0.74417
Time: 3.6 Velocity: 0.832378 Position: -0.702551
Time: 3.65 Velocity: 1.09311 Position: -0.647895
Time: 3.7 Velocity: 1.33514 Position: -0.581138
Time: 3.75 Velocity: 1.5538 Position: -0.503448
Time: 3.8 Velocity: 1.74461 Position: -0.416217
Time: 3.85 Velocity: 1.90345 Position: -0.321045
Time: 3.9 Velocity: 2.02668 Position: -0.219711
Time: 3.95 Velocity: 2.11141 Position: -0.11414
Time: 4 Velocity: 2.15555 Position: -0.00636285
Time: 4.05 Velocity: 2.15801 Position: 0.101538
Time: 4.1 Velocity: 2.11874 Position: 0.207475
Time: 4.15 Velocity: 2.0387 Position: 0.309409
Time: 4.2 Velocity: 1.91985 Position: 0.405402
Time: 4.25 Velocity: 1.76503 Position: 0.493653
Time: 4.3 Velocity: 1.57778 Position: 0.572542
Time: 4.35 Velocity: 1.36216 Position: 0.640651
Time: 4.4 Velocity: 1.12264 Position: 0.696783
Time: 4.45 Velocity: 0.863863 Position: 0.739976
Time: 4.5 Velocity: 0.590567 Position: 0.769504
Time: 4.55 Velocity: 0.307511 Position: 0.78488
Time: 4.6 Velocity: 0.0194268 Position: 0.785851
Time: 4.65 Velocity: -0.268974 Position: 0.772402
Time: 4.7 Velocity: -0.552981 Position: 0.744753
Time: 4.75 Velocity: -0.827865 Position: 0.70336
Time: 4.8 Velocity: -1.08887 Position: 0.648916
Time: 4.85 Velocity: -1.33126 Position: 0.582354
Time: 4.9 Velocity: -1.55035 Position: 0.504836
Time: 4.95 Velocity: -1.74166 Position: 0.417753
Time: 5 Velocity: -1.90106 Position: 0.3227
Press any key to continue . . .

There you go. All the time units are seconds velocity is m/s and position is radians from equilibrium, positive to the right, negative to the left.


Yeah there are several things i can do, like you mentioned Integral in order to streamline and increase the speed of the code, but at this point I'm just happy that its working. I'll probably spend some time this weekend streamlining it as best i can. However i prefer to leave g as an input value, only because then it has the added (although possibly pointless) functionality of comparing pendulum behavior under different amounts of gravity. The computing g/l once though is a good idea, though it probably won't save much time (10 seconds at 0.01 time intevals only takes about 3 seconds to compute).

And agian, thanks for the help.

 

[Integral]

When I was actively doing this sort of calculation,it was on a 8 bit 1Mhz Apple II computer.(yes, it was 20yrs ago!) There were several models that I han to let run over a weekend to complete.

While the modern computers are incredibly capable I would like to think that leaner more efficient codes are simply more pleasing to create.

You might note the deviation of the period from that predicted by the usual equation. Try running your program with a smaller starting angle to see how the period varies with initial position.

I have been looking for some error analysis of this method, have not yet found it, I believe that the error goes like h^4. One way to get a feel for the errors is to complete the same run with a range of step sizes, then observe which digits change with the change in step.

 

[nbo10]

What method are you looking for error analysis? Runge-Kutta?? second order RK is h^4, fourth order RK is h^8. but you do a pay price. The code I posted was just a outline of a way what you need to do. If you run the that I posted you can see that something isn't quite right?

JMD

 

[Integral]

Why do you think there is a problem?

Franz,
Could you run your program with an initial position of .1 and post the results? TIA

 

[franznietzsche]

Originally posted by Integral
Why do you think there is a problem?

Franz,
Could you run your program with an initial position of .1 and post the results? TIA

Sure thing:

Enter a value for g: -9.8

Enter an initial angular position:0.1

Enter the pendulum length: 1

Enter the time interval between predictions: 0.05

Enter run-time length of simulation: 5

Press any key to continue . . .
Time: 0 Velocity: 0 Position: 0.1
Time: 0.05 Velocity: -0.0386838 Position: 0.0980658
Time: 0.1 Velocity: -0.0766206 Position: 0.0942348
Time: 0.15 Velocity: -0.113078 Position: 0.0885809
Time: 0.2 Velocity: -0.14735 Position: 0.0812134
Time: 0.25 Velocity: -0.178775 Position: 0.0722746
Time: 0.3 Velocity: -0.206745 Position: 0.0619374
Time: 0.35 Velocity: -0.230717 Position: 0.0504016
Time: 0.4 Velocity: -0.250226 Position: 0.0378903
Time: 0.45 Velocity: -0.264894 Position: 0.0246456
Time: 0.5 Velocity: -0.274435 Position: 0.0109238
Time: 0.55 Velocity: -0.278664 Position: -0.00300937
Time: 0.6 Velocity: -0.277499 Position: -0.0168843
Time: 0.65 Velocity: -0.270962 Position: -0.0304324
Time: 0.7 Velocity: -0.259181 Position: -0.0433915
Time: 0.75 Velocity: -0.242384 Position: -0.0555107
Time: 0.8 Velocity: -0.220898 Position: -0.0665556
Time: 0.85 Velocity: -0.19514 Position: -0.0763127
Time: 0.9 Velocity: -0.16561 Position: -0.0845931
Time: 0.95 Velocity: -0.132878 Position: -0.091237
Time: 1 Velocity: -0.0975793 Position: -0.096116
Time: 1.05 Velocity: -0.0603957 Position: -0.0991358
Time: 1.1 Velocity: -0.0220456 Position: -0.100238
Time: 1.15 Velocity: 0.0167301 Position: -0.0994016
Time: 1.2 Velocity: 0.0551828 Position: -0.0966424
Time: 1.25 Velocity: 0.0925698 Position: -0.0920139
Time: 1.3 Velocity: 0.128169 Position: -0.0856055
Time: 1.35 Velocity: 0.161292 Position: -0.0775409
Time: 1.4 Velocity: 0.191297 Position: -0.0679761
Time: 1.45 Velocity: 0.217604 Position: -0.0570958
Time: 1.5 Velocity: 0.239704 Position: -0.0451107
Time: 1.55 Velocity: 0.257166 Position: -0.0322524
Time: 1.6 Velocity: 0.269651 Position: -0.0187698
Time: 1.65 Velocity: 0.276918 Position: -0.00492394
Time: 1.7 Velocity: 0.278824 Position: 0.00901727
Time: 1.75 Velocity: 0.275333 Position: 0.0227839
Time: 1.8 Velocity: 0.266513 Position: 0.0361095
Time: 1.85 Velocity: 0.252534 Position: 0.0487362
Time: 1.9 Velocity: 0.233669 Position: 0.0604197
Time: 1.95 Velocity: 0.210284 Position: 0.0709339
Time: 2 Velocity: 0.182833 Position: 0.0800756
Time: 2.05 Velocity: 0.151848 Position: 0.087668
Time: 2.1 Velocity: 0.117928 Position: 0.0935644
Time: 2.15 Velocity: 0.08173 Position: 0.0976509
Time: 2.2 Velocity: 0.0439535 Position: 0.0998485
Time: 2.25 Velocity: 0.00532819 Position: 0.100115
Time: 2.3 Velocity: -0.0334 Position: 0.0984449
Time: 2.35 Velocity: -0.0714832 Position: 0.0948708
Time: 2.4 Velocity: -0.108186 Position: 0.0894615
Time: 2.45 Velocity: -0.142799 Position: 0.0823216
Time: 2.5 Velocity: -0.174652 Position: 0.0735889
Time: 2.55 Velocity: -0.20313 Position: 0.0634324
Time: 2.6 Velocity: -0.22768 Position: 0.0520484
Time: 2.65 Velocity: -0.247826 Position: 0.0396571
Time: 2.7 Velocity: -0.263178 Position: 0.0264982
Time: 2.75 Velocity: -0.273436 Position: 0.0128264
Time: 2.8 Velocity: -0.278402 Position: -0.0010937
Time: 2.85 Velocity: -0.277979 Position: -0.0149926
Time: 2.9 Velocity: -0.272174 Position: -0.0286014
Time: 2.95 Velocity: -0.261102 Position: -0.0416564
Time: 3 Velocity: -0.244977 Position: -0.0539053
Time: 3.05 Velocity: -0.224112 Position: -0.0651109
Time: 3.1 Velocity: -0.198913 Position: -0.0750565
Time: 3.15 Velocity: -0.169867 Position: -0.0835499
Time: 3.2 Velocity: -0.137539 Position: -0.0904268
Time: 3.25 Velocity: -0.102553 Position: -0.0955545
Time: 3.3 Velocity: -0.0655866 Position: -0.0988338
Time: 3.35 Velocity: -0.0273532 Position: -0.100201
Time: 3.4 Velocity: 0.0114084 Position: -0.099631
Time: 3.45 Velocity: 0.0499498 Position: -0.0971336
Time: 3.5 Velocity: 0.0875265 Position: -0.0927572
Time: 3.55 Velocity: 0.123413 Position: -0.0865866
Time: 3.6 Velocity: 0.156914 Position: -0.0787409
Time: 3.65 Velocity: 0.187384 Position: -0.0693717
Time: 3.7 Velocity: 0.214231 Position: -0.0586601
Time: 3.75 Velocity: 0.236935 Position: -0.0468134
Time: 3.8 Velocity: 0.255056 Position: -0.0340606
Time: 3.85 Velocity: 0.268242 Position: -0.0206485
Time: 3.9 Velocity: 0.276235 Position: -0.00683671
Time: 3.95 Velocity: 0.278882 Position: 0.00710741
Time: 4 Velocity: 0.276131 Position: 0.0209139
Time: 4.05 Velocity: 0.268034 Position: 0.0343156
Time: 4.1 Velocity: 0.25475 Position: 0.0470531
Time: 4.15 Velocity: 0.236536 Position: 0.0588799
Time: 4.2 Velocity: 0.213747 Position: 0.0695673
Time: 4.25 Velocity: 0.186824 Position: 0.0789085
Time: 4.3 Velocity: 0.15629 Position: 0.086723
Time: 4.35 Velocity: 0.122735 Position: 0.0928598
Time: 4.4 Velocity: 0.0868096 Position: 0.0972002
Time: 4.45 Velocity: 0.0492072 Position: 0.0996606
Time: 4.5 Velocity: 0.0106544 Position: 0.100193
Time: 4.55 Velocity: -0.028104 Position: 0.0987881
Time: 4.6 Velocity: -0.0663198 Position: 0.0954721
Time: 4.65 Velocity: -0.103255 Position: 0.0903094
Time: 4.7 Velocity: -0.138195 Position: 0.0833996
Time: 4.75 Velocity: -0.170465 Position: 0.0748764
Time: 4.8 Velocity: -0.199441 Position: 0.0649043
Time: 4.85 Velocity: -0.22456 Position: 0.0536763
Time: 4.9 Velocity: -0.245336 Position: 0.0414095
Time: 4.95 Velocity: -0.261366 Position: 0.0283412
Time: 5 Velocity: -0.272338 Position: 0.0147243
Press any key to continue . . .

There is some visible error by the first half period from what would be expected, but it is barely measurable in terms of the kinds of apparatus i have access to (meaning a digital video camera with frame by frame play back). Also the period for a pendulum of this length should be about 2.08 seconds but at 1.1 seconds it still had not quite changed directions.

Also i agree about the creating leaner code being a generally better idea, even if just for the sake of practice. I'll probably declare all the pkone-four and vkone-four functions to be inline reducing the number of function calls per iteration to two instead of the muhc higher number it is now.

 

[Integral]

Franz,

There must remain a significant error in your code. The last run should match the small angle approximation solution of

$$\theta (t)= .1 cos( \frac {gt} l)$$

This says the max velocity should .1 g/l. if you look at the run I did in Excel with a simple Euler method you can verify that. Also you should have a period of 2.007 seconds, this corresponds to my Euler run but not your last run.

Unfortunatly there remains a problem in your code. I cannot see it yet but will return to this shortly

(My wifes car has a flat, seems that she thinks that problem has priority! The Nerve!)

 

[franznietzsche]

Originally posted by Integral
Franz,

There must remain a significant error in your code. The last run should match the small angle approximation solution of

$$\theta (t)= .1 cos( \frac {gt} l)$$

This says the max velocity should .1 g/l. if you look at the run I did in Excel with a simple Euler method you can verify that. Also you should have a period of 2.007 seconds, this corresponds to my Euler run but not your last run.

Unfortunatly there remains a problem in your code. I cannot see it yet but will return to this shortly

(My wifes car has a flat, seems that she thinks that problem has priority! The Nerve!)

Using conservation of energy and a starting position of 0.1 radians i calculate a maximum veolcity of 0.312919 m/s. then running the program again i get:

Enter a value for g: -9.8

Enter an initial angular position:0.1

Enter the pendulum length: 1

Enter the time interval between predictions: .100354

Enter run-time length of simulation: 2.00709

Press any key to continue . . .
Time: 0 Velocity: 0 Position: 0.1
Time: 0.100354 Velocity: -0.0619174 Position: 0.0937863
Time: 0.200708 Velocity: -0.11999 Position: 0.0817448
Time: 0.301062 Velocity: -0.170611 Position: 0.0646233
Time: 0.401416 Velocity: -0.210634 Position: 0.0434854
Time: 0.50177 Velocity: -0.237567 Position: 0.0196446
Time: 0.602124 Velocity: -0.249735 Position: -0.00541737
Time: 0.702478 Velocity: -0.24638 Position: -0.0301426
Time: 0.802832 Velocity: -0.22771 Position: -0.0529941
Time: 0.903186 Velocity: -0.194888 Position: -0.0725519
Time: 1.00354 Velocity: -0.149957 Position: -0.0876006
Time: 1.10389 Velocity: -0.0957117 Position: -0.0972057
Time: 1.20425 Velocity: -0.0355232 Position: -0.100771
Time: 1.3046 Velocity: 0.0268709 Position: -0.098074
Time: 1.40496 Velocity: 0.0875967 Position: -0.0892833
Time: 1.50531 Velocity: 0.142883 Position: -0.0749444
Time: 1.60566 Velocity: 0.189295 Position: -0.0559479
Time: 1.70602 Velocity: 0.223946 Position: -0.0334741
Time: 1.80637 Velocity: 0.244679 Position: -0.00891951
Time: 1.90673 Velocity: 0.250204 Position: 0.0161895
Time: 2.00708 Velocity: 0.240176 Position: 0.0402922
Press any key to continue . . .

the max velocity should occur at 0.50177 sec, but it doesn't. the question then is whether the eroor is from the code, or the Runge-Kutta algorithm itself (since the error would be squared due to an application of the algorithm to its own results.)

Edit: Pff, everyuone knows Newtonian mechanics are WAY more important than auto mechanics. silly her..hehe.

 

[franznietzsche]

Also i just noticed something else in the errors:
Enter a value for g: -9.8

Enter an initial angular position:0.785398

Enter the pendulum length: 1

Enter the time interval between predictions: 0.0100354

Enter run-time length of simulation: 2.50886

Press any key to continue . . .



Time: 1.53542 Velocity: 2.30635 Position: -0.129727
Time: 1.54545 Velocity: 2.31847 Position: -0.10646
Time: 1.55549 Velocity: 2.32843 Position: -0.0830934
Time: 1.56552 Velocity: 2.3362 Position: -0.0596486
Time: 1.57556 Velocity: 2.34179 Position: -0.0361479
Time: 1.58559 Velocity: 2.34517 Position: -0.0126131
Time: 1.59563 Velocity: 2.34635 Position: 0.0109335
Time: 1.60566 Velocity: 2.34533 Position: 0.0344698
Time: 1.6157 Velocity: 2.3421 Position: 0.0579737

Time: 2.00708 Velocity: 0.795103 Position: 0.74091
Time: 2.01712 Velocity: 0.731119 Position: 0.748247
Time: 2.02715 Velocity: 0.666608 Position: 0.754937
Time: 2.03719 Velocity: 0.60162 Position: 0.760974
Time: 2.04722 Velocity: 0.536204 Position: 0.766355
Time: 2.05726 Velocity: 0.470407 Position: 0.771076
Time: 2.06729 Velocity: 0.404279 Position: 0.775133
Time: 2.07733 Velocity: 0.337867 Position: 0.778524
Time: 2.08736 Velocity: 0.271217 Position: 0.781246
Time: 2.0974 Velocity: 0.204378 Position: 0.783297
Time: 2.10743 Velocity: 0.137396 Position: 0.784676
Time: 2.11747 Velocity: 0.0703191 Position: 0.785381
Time: 2.1275 Velocity: 0.00319285 Position: 0.785413
Time: 2.13754 Velocity: -0.0639356 Position: 0.784772
Time: 2.14758 Velocity: -0.13102 Position: 0.783457
Time: 2.15761 Velocity: -0.198012 Position: 0.78147

the position should max out at 2.00708 seconds, but it doesn't, it keeps rising until 2.1275 seconds, at which point it is only .001 radians off from what the maximum position should be. Also the highest velocity listed is 2.34635 m/s, which is only 0.04963 m/s less than the predicted maximum veolcity. The inaccuracy exists in that these are predicted at the wrong times, but the actual predicted maximum and minimum values of position and velocity are very close to what they should be.

 

[franznietzsche]

Also another thought just occered to me. the Newtonian equation for the period of a pendulum $$T = 2*pi*sqrt(l/g)$$ is only an approximation and only accurate for small displacements. Now i noticed that with the 0.1 starting position the program predicted a shorter period that the Newtonian equation, but for the 0.785398 radian initial displacement it predicted a longer period.

Now is it possible that the computer program is following the actual period more closely than Newton's equation? Does anyone here know what the real equation for the period of a pendulum (or other harmonic oscillator) is?

 

[Integral]

That is why I wanted to see the run with initial position of .1, it should be modeled very nicly by the closed form approximate solution. For that run you would expect the period to be spot on. That is not true for the ~.7 initial postion, we would expect to see some variation at that amplitude. If we cannot get good results with .1 amplitude, then we cannot get meaningful results anywhere.


edit: yes, English is my native language, believe it or not!