Predicting Pendulum Position: Differential Equations in C++ Code

[Integral]

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!)

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

 

[franznietzsche]

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.

 

[franznietzsche]

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

Enter run-time length of simulation: 2.00709

Press any key to continue . . .

Time: 0.481699 Velocity: -0.303909 Position: 0.00845068
Time: 0.491735 Velocity: -0.304701 Position: 0.00539289
Time: 0.50177 Velocity: -0.305206 Position: 0.00233003
Time: 0.511805 Velocity: -0.305424 Position: -0.000735027
Time: 0.521841 Velocity: -0.305355 Position: -0.00379939

Press any key to continue . . .

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.

 

[franznietzsche]

i wrote another program to calculate the period using the previous program's Runge-Kutta algoithm. it is:

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

 

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

 

[franznietzsche]

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.

 

[Integral]

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!

k1 = h * Velocity_old
l1 = h * RHS(L, Position_old)
k2 = h * (Velocity_old + (l1 / 2))
l2 = h * RHS(L, Position_old + (k1 / 2))
k3 = h * (Velocity_old + (l2 / 2))
l3 = h * RHS(L, Position_old + k2 / 2)
k4 = h * (Velocity_old + l3)
l4 = h * RHS(L, Position_old + k4)

      Velocity_new = Velocity_old + (l1 + 2 * l2 + 2 * l3 + l4) / 6
      Position_new = Position_old + (k1 + 2 * k2 + 2 * k3 + k4) / 6

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!

 

[franznietzsche]

Originally posted by Integral

In this code snipit RHS is for Right Hand Side, it is your Acc.

QUOTE]

??

 

[Integral]

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?

 

[nbo10]

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.

JMD

 

[Integral]

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.

 

[nbo10]

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.

JMD

 

[nbo10]

Euler-Kromer method

JMD

 

[Integral]

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.

 

[Integral]

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.

 

[nbo10]

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.

t[i+1] = t[i] + dt;
omega[i+1] = omega[i] - (g/l) * theta[i] * dt;
theta[i+1] = theta[i] + omega[i] * dt;
Energy[i+1] = g*l*(1 - Math.Cos(3.1415*theta[i])/180 ) + 0.5 * omega[i]*omega[i];

This doesn't work, Euler method.

t[i+1] = t[i] + dt;
omega[i+1] = omega[i] - (g/l) * theta[i] * dt;
theta[i+1] = theta[i] + omega[i+1] * dt;
Energy[i+1] = g*l*(1 - Math.Cos(3.1415*theta[i])/180 ) + 0.5 * omega[i]*omega[i];

This does work, Euler-Kromers method.

JMD

 

[Integral]

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.

 

[nbo10]

Using the closed form of the solution enhances the problem of energy conservation.

do {
i++;
t[i+1] = t[i] + dt;
omega[i+1] = omega[i] - (g/l) * Math.Sin(3.1415*theta[i]/180)* dt;
theta[i+1] = theta[i] + omega[i] * dt;
Energy[i+1] = g*l*(1 - Math.Cos(3.1415*theta[i])/180 ) + 0.5 *l*l* omega[i]*omega[i];

omega1[i+1] = omega1[i] - (g/l) * Math.Sin(3.1415*theta1[i]/180) * dt;
theta1[i+1] = theta1[i] + omega1[i+1] * dt;
Energy1[i+1] = g*l*(1 - Math.Cos(3.1415*theta1[i])/180 ) + 0.5 * l*l*omega1[i]*omega1[i];

while (t[i+1] <= 10*per);
thetaplot.PlotXY(t,Energy1);
omegaplot.PlotXY(t,Energy);

theta, omega, energy, use Euler method.
theta1, omega1, energy1 use Euler-Kromer method.
Here is the result a plot of energy vs time.

I believe Numerical Analysis, 6th Ed. Talks about the problems of Eulers method. I can't be for sure because, I don't have my copy handy.

JMD

 

[Integral]

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.

 

[nbo10]

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.

JMD

 

[Integral]

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.

 

[Integral]

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.

 

[Stingray]

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

 

[Integral]

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? :)

 

[krab]

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.

 

[Integral]

Nice link Krab. What I fail to see there is any formal error analysis of the methods. The claim that the "failure to conserve energy" property of Euler is characteristic to periodic functions is not proven. They only show 1 step size, what would happen if they choose a smaller step? What would happen if they chose a DIFFERENT function. The local discretization error of Euler is ~step size, this holds for ALL applications of Euler,not just periodic.

Here is an expression derived in Conte and de Boor which relates step size and local error.

$$|f_y(x,y)|\leq L$$
$$|y^{\prime\prime}(x)|<Y$$
$$|{error}_n|= \frac{hY}{2L}(e^{(x_n-x_0)L}-1)$$

He points out that this is a "Upper bound rather then a realistic bound" so you may be able to achieve the same error with a smaller step. But a step size computed with this formula will yield the desired error at the specified end point.

Notice that this error is indeed monotonic. So for the pendulum it shows up as a steadily increasing amplitude. Some one with no deeper vision could indeed say that it looks like energy is not conserved. But that is failing to look at the mechanism of the tool being applied. (Failure: Car won't run, solution: Replace engine, root cause: out of gas) Further, if you were apply this method to a different system (perhaps financial) would it still fail to conserve energy? Of course not, because energy is not involved, but it would still show a steady increase in the final solution. So to say Euler fails to conserve energy is to broad a statement. It can only be make for specific models. While when I say Euler has local discretization of ~step this applies to ALL applications of Euler, and is the root cause of the failure to conserve energy phenomena in the pendulum problem.

 

[franznietzsche]

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

I just did two separate runs one using Euler and one using Euler-Cromer both with dt = 0.01. Euler-Cromer was very accurate at 100 seconds and gave the correct time for the equilibrium points to within a factor of 0.03 of the period predicted by Newton's equation. When i did the run with Euler's Method however soemthing interesting happened. The period became longer, however energy was conserved to the same degree it had been with Euler-Cromer. I noted no difference in amplitude of either velocity or position, though the period was longer, but only marginally so. At 100 seconds Eulers method was about a second behind the Euler-Cromer method.

 

[Integral]

I have seen similer results, where the amplitude (energy) remained consitent but the period was incorrect. It did not seem to creep out, it was just wrong. Not real sure why that happens. I wondered if I was writing to the wrong data file.