Predicting Pendulum Position: Differential Equations in C++ Code

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