The discussion centers on a C++ program designed to predict the position of a pendulum using differential equations, specifically the motion equations. The user experiences issues with the simulation, where the pendulum's amplitude and period increase unrealistically, leading to nonsensical results such as a position of -18 radians after 10 seconds. Contributors point out that the original equations used were incorrect, suggesting the need for the correct form of angular acceleration, which should include a negative sign. The conversation also highlights the importance of using numerical methods like the Runge-Kutta algorithm for better accuracy in approximating the pendulum's motion. Ultimately, the user is encouraged to refine their implementation to achieve more realistic oscillatory behavior.
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!)
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
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.
Last edited:
#33
franznietzsche
1,503
6
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
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.
#34
franznietzsche
1,503
6
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?
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!
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!)
\
At .1 initial postion the difference between this equation and our RK model should be on the order of h^4, we should agree to 3-4 decimals.
BTW g in this model must be positive, its sign has been factored in already.
Last edited:
#37
franznietzsche
1,503
6
Originally posted by Integral \
At .1 initial postion the difference between this equation and our RK model should be on the order of h^4, we should agree to 3-4 decimals.
BTW g in this model must be positive, its sign has been factored in already.
the sign isn't built-in in my runge-kutta algorith though because the angular is acceleration is (g/l)*sin(theta). SO i input it as a negative.
Also I'm going to run the program for many variations of g, l, and initial position to try and see if the period that the program gives is a consistent function of those values and how close it is to Newton's equation. It is certainly more likely that the model is simply flawed, but it is worthwhile to examine it more closely, because it is also known that Newton's equation is overly simplistic, and because the period was smaller in the model at a low amplitude, and larger in the model for a larger amplitude i think it is possible even if unlikely that the model is following what the real equation predicts within the predicted error.
#38
franznietzsche
1,503
6
ran a new simulation for h=0.0100354 and initial position 0.1:
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.0100354
the actual calculation for the equilibrium point is t = 0.5017724078863 seconds, and the max velocity would be -0.3129. This h value is ten times smaller than the previous one and much more accurate. I think now that the problem is that error is doubled, not squared or some other change when using the algorithm on itself.
#39
franznietzsche
1,503
6
i wrote another program to calculate the period using the previous program's Runge-Kutta algoithm. it is:
Code:
#include<iostream.h>
#include<math.h>
#include<stdlib.h>
#include<pendulumfunctions.h>
int main()
{
g = -9.8;
cout << "Enter the length of the pendulum: ";
cin >> l;
gl = g/l;
cout << "Enter the initial angular position: ";
cin >> ang;
cout << "Enter duration of simulation: ";
cin >> z;
iang = ang;
y = ang;
w = 0;
av = 0;
t = 0;
h=0.00001;
while (t<z)
{
ang = y;
av = w;
velocity(gl, ang);
position(gl, ang);
if ((iang-ang)<0.000000001)
{
cout << "Time: " << t << "\tVelocity: " << velocity(gl, ang) << "\tPosition: " << position(gl, ang) << endl;
system("PAUSE");
}
t = t + h;
}
system("PAUSE");
return 0;
}
[B]Pendulumfunctions.h:[/B]
#include<math.h>
typedef double dou;
dou g, l, ang, iang, av, aa, t, h, w, y, v, z, gl;
dou aac(dou gl, dou ang) //Function for angular acceleration, radians/sec
{
dou aa;
aa = gl*sin(ang);
return (aa);
}
inline dou vkone()
{
return(aac(gl, ang));
}
inline dou vktwo()
{
return(aac(gl, (ang + (h/2)*vkone())));
}
inline dou vkthree()
{
return(aac(gl, (ang + (h/2)*vktwo())));
}
inline dou vkfour()
{
return(aac(gl, (ang + h*vkthree())));
}
dou ave(dou gl, dou ang) //Uses Runge-Kutta Algorithm for angular velocity
{
if (t==0)
return av;
else
{
w = av +(h/6)*(vkone() + 2*vktwo() + 2*vkthree() + vkfour());
return w;
}
}
dou velocity(dou gl, dou ang)
{
v = ave(gl, ang)*l;
return (v);
}
inline dou pkone()
{
return(ave(gl, ang));
}
inline dou pktwo()
{
return(ave(gl, ang + (h/2)*pkone()));
}
inline dou pkthree()
{
return(ave(gl, ang + (h/2)*pktwo()));
}
inline dou pkfour()
{
return(ave(gl, (ang + h*pkthree())));
}
dou position(dou gl, dou ang) //Uses Runge Kutta algorithm for position
{
if (t==0)
return iang;
else
{
y = ang + (h/6)*(pkone() + 2*pktwo() + 2*pkthree() + pkfour());
return y;
}
}
The header file pendulumfunctions.h is unchanged between the two programs, its the main file that differs. here are the period results for 0.1 radians starting position and 0.785398 starting position:
Enter the length of the pendulum: 1
Enter the initial angular position: 0.1
Enter duration of simulation: 8 Time: 0 Velocity: 0 Position: 0.1
Press any key to continue . . .
Time: 1e-05 Velocity: -9.7832e-06 Position: 0.1
Press any key to continue . . .
Time: 2e-05 Velocity: -1.95664e-05 Position: 0.1
Press any key to continue . . .
Time: 3e-05 Velocity: -2.93496e-05 Position: 0.1
Press any key to continue . . .
Time: 4e-05 Velocity: -3.91328e-05 Position: 0.1
Press any key to continue . . .
Time: 5e-05 Velocity: -4.8916e-05 Position: 0.1
Press any key to continue . . .
Time: 2.00836 Velocity: 3.33546e-05 Position: 0.1
Press any key to continue . . .
Time: 2.00837 Velocity: 2.35714e-05 Position: 0.1
Press any key to continue . . .
Time: 2.00838 Velocity: 1.37882e-05 Position: 0.1
Press any key to continue . . . Time: 2.00839 Velocity: 4.00496e-06 Position: 0.1
Press any key to continue . . .
Time: 2.0084 Velocity: -5.77824e-06 Position: 0.1
Press any key to continue . . .
Time: 2.00841 Velocity: -1.55614e-05 Position: 0.1
Press any key to continue . . .
Time: 2.00842 Velocity: -2.53446e-05 Position: 0.1
Press any key to continue . . .
Time: 2.00843 Velocity: -3.51278e-05 Position: 0.1
Press any key to continue . . .
Time: 2.00844 Velocity: -4.4911e-05 Position: 0.1
Press any key to continue . . .
and for 0.785398 (pi/4 or 45 degrees):
Enter the length of the pendulum: 1
Enter the initial angular position: 0.785398
Enter duration of simulation: 8 Time: 0 Velocity: 0 Position: 0.785398
Press any key to continue . . .
Time: 1e-05 Velocity: -6.92941e-05 Position: 0.785398
Press any key to continue . . .
Time: 2e-05 Velocity: -0.000138588 Position: 0.785398
Press any key to continue . . .
Time: 2.08735 Velocity: 6.80124e-05 Position: 0.785398
Press any key to continue . . . Time: 2.08736 Velocity: -1.28162e-06 Position: 0.785398
Press any key to continue . . .
Time: 2.08737 Velocity: -7.05757e-05 Position: 0.785398
Press any key to continue . . .
Time: 2.08738 Velocity: -0.00013987 Position: 0.785398
Press any key to continue . . .
I bolded the times that gave the lowest velocities. Now these two use a very low h value for the algorithm, h = 0.00001 which causes the period predicted on the first one to be within .001 of the predicted period. The period for the larger amplitude is farther from the predicted, but this is too be expected as Newton's equation becomes invalid for large initial amplitudes.
Based upon this data it looks like the inaccuracy noted earlier was caused not by any flaws in the algorithm (the two programs use the exact same algorithm) but by using an h value that was altogether too large. So the predictions made by the program that gives the position are accurrate given a low enough value of h.
http://home.comcast.net/~rossgr1/Math/Eulererror.pdf is a PDF showing the deviation of a simple Euler method from the closed form solution. With a similer step size the RK method should be a better approximation, it is not. I have coded in VB both a Euler method (it produced the numbers in this pdf) and an RK method. My RK is produceing similar results as yours. I am staring at the equations, and think I am close to finding the problem in our implementation.
Last edited by a moderator:
#41
franznietzsche
1,503
6
Originally posted by Integral http://home.comcast.net/~rossgr1/Math/Eulererror.pdf is a PDF showing the deviation of a simple Euler method from the closed form solution. With a similer step size the RK method should be a better approximation, it is not. I have coded in VB both a Euler method (it produced the numbers in this pdf) and an RK method. My RK is produceing similar results as yours. I am staring at the equations, and think I am close to finding the problem in our implementation.
I think the problem is the application of the RK algorithm to its ownresults. I've been doing regression analysis for the past few hours of the predicted periods, and for small enough theta they match Newton's equation near exactly. What i have is a quadratic (with theta being the variable that is squared) with very small (about .06159L+0.078092 and 0.019259L+0.017962 respectively) a and b coefficients, and the c value is Newton's equation for the period.
Ok, http://home.comcast.net/~rossgr1/Math/Rkerror.PDF is the way an R-K method SHOULD look! Note that this run has the same step size as the Euler method I posted earlier, but has at least 2 more good digits!
I finally found the section in Conte and Deboor dealing with R-K and SYSTEMS of equations. You cannot compute the corrections seperatly, notice in the above code the ls are corrections for velocity but they depend on the ks and visa versa.
In this code snipit RHS is for Right Hand Side, it is your Acc.
This works.
Whew... Being beating my head on this all day!
Last edited by a moderator:
#43
franznietzsche
1,503
6
Originally posted by Integral
In this code snipit RHS is for Right Hand Side, it is your Acc.
Function RHS(Length, Position)
g = 9.80665
RHS = -g * Sin(Position) / Length
End Function
Is that better? :)
(I carry the negitive with the function.)
Opps... Make that it is your aac
Edit one more time!:
I see that you are factoring your step out of the k2 (your Vkone ) . Can't do that! note that since k2 depends on l1 if you do not have the factor of h in l1 , k2 will be wrong.
That make sense?
Last edited:
#45
nbo10
416
5
I've finally had some time to think about this. The euler method doesn't conserve energy. You can see this by graphing theta vs t. you'll see theta_max increases with each period. YOu can use Euler-Cromer method which is the next easiest method.
nbo10,
The Euler method works fine. I do not see any increase in amplitude with time. As well as the fact that it has been in use for years and has always worked fine, considering its inherent errors. Perhaps you need to examine your implementation.
#47
nbo10
416
5
In general Eulers methods fine, but not when you use it to solve the pendulum problem. Calculate the energy and graph the results. The attached gif is what you'll see. YOu can decrease the error in Eulers method by decreaseing dt, but you still won't conserve energy.
Once again, it is your implementation that is at fault. Not Euler.
Since the Euler method is simply a discretization of a derivative, it cannot effect the energy of the system. However, careless implementation can result in garbage, such as you are seeing.
http://home.comcast.net/~rossgr1/Math/Eulerc.PDF is a chart of my implementation of the Euler method. There is NO tendency to increase in amplitude. Please stop blaming the method for your implementation errors.
Last edited by a moderator:
#51
nbo10
416
5
I'm not going to sit here and argue about this, I have too much to do. For periodic functions Eulers method is unstable. Run YOUR program using Eulers method for 100 periods check the results. It's unstable. THis is my last post about this subject.
Please provide some references concerning the instability of Euler for periodic functions.
First of note that you are using the small angle approximation to the pendulum. It is not necessary to apply numerical methods to that problem. It has a nice closed form solution which I have already posted in this thread.
I am not sure what your problem is. Euler could well fail after 100 periods simply due to accumulated error. It is prone to that. R-K is much better in that respect. To the best of my knowledge, Euler is not any more unstable for Periodic functions then it is for any other class of function.
If you have some hard core references please provide them. If not then please speak carefully when expressing such opinions.
#53
nbo10
416
5
Using the closed form of the solution enhances the problem of energy conservation.
I can find no reference to Kromer in my books. But looking at what you are doing it is clear that your Kromers method is simply the correct implementation of Euler's method to a SYSTEM of equations. Your first code (what you call simply Euler's) is not the way to solve a system. While the second is. So yes your first program will not work.
Edit:
There is no doubt that there are inherent troubles with Euler, even if implemented correctly. But it is no more unstable or inaccurate for Periodic solutions as any other 2nd order equation.
Last edited:
#55
nbo10
416
5
It's Cromer and not Kromer, my mistake. Every book I've read has Euler as the way that I have it above.
You can show that in Euler-Cromer method (for oscillatory motion over a period) the Energy increase is ~dt^3, while Euler method energy is increase is ~dt.
Remember that this is a system of equations. Look at the variation of the Runga Kutta method that was need to correctly model a system of equations. The same follows for Euler's. The trouble comes with how to best solve the system as much as the method used.
Using the closed form of the solution enhances the problem of energy conservation.
Could you please expand on this?
I cannot see how you can meaingfully relate discretiontion error to energy conservation. If there were something physically incorrect, as you imply, then it would not matter how small a step size you take, the method would still "not conserve energy". The fact is Euler has error ~h (or dt as you call it) So what you say is trivial and has nothing to do with energy or periodic functions.
Integral, numerical analysis involves a lot more than using Taylor series to show roundoff error.
Certain methods with apparently identical roundoff can be better or worse at solving specific problems. Periodic motion tends to rack up large numerical errors after many periods, so naive implementations won't last long. nbo's labeling of "Euler" and "Euler-Cromer" methods is standard, and it is easy to show that Euler's method leads to monotonically increasing energy (if I remember right - I don't feel like writing anything right now), whereas Euler-Cromer does not.
Codes for more difficult systems like hydrodynamics often have conserved quantities explicitly built into the integration scheme. Things are much more stable that way. Runge-Kutta is useful for lots of things, but its not always the best tool.
Anyways, the pendulum can be solved analytically without the small angle approximation in terms of elliptic integrals. The properties of these functions are commonly tabulated in books, you can plot them in mathematica, etc (just like trig functions...).
My posts are based on the error analysis given in Elementary Numerical Analysis by Conte and de Boor and A Survey of Numerical Methods by Young and Gregory.
Certainly a compounding discretization error will lead to diverging amplitudes. You can call this failure to conserve energy if you wish, I will call it discretization error. Mainly because it can be controlled by choosing a smaller step size. Though why anybody would want to use Euler's or slight variations is not clear. Certainly there are better methods then R-K but it works.
I am surprised that the Cromer variation is not named in any of my texts, though comparing to Nob code, the Euler's numbers I generated were with the Cromers variation. I did it that way because it felt good!
One thing I do know about Numerical Methods is that you must watch your numbers like a hawk. The line between a stable solution and garbage is very thin. That is the reason I wanted to see some small amplitude numbers, those can be verified independently. If a model can be shown to be good for one set of IC it will also be good for other IC which are not hugely different.
I suppose you could solve the pendulum with elliptic integrals... What fun would that be? :)
Look here
for a comparison of Euler, Euler-Cromer, Runge-Kutta methods. It includes Matlab scripts.
In particular, notice that the straight Euler method gives energy error that grows exponentially; so whoever claimed first that Euler method is unstable for periodic motion appears to be correct.
I've done a lot of this kind of problem and have always used RK; it has never let me down.
