Pendulum simulation problem

george7378

DON'T PANIC
Addon Developer
Donator
Joined
Jun 26, 2009
Messages
1,045
Reaction score
0
Points
36
Hi everyone,

I made a small program to simulate a pendulum last night and it seems to work OK apart from one thing: the pendulum gains momentum as it swings. First off, here's the code I'm using:

Code:
accAng = -g*sin(posAng)/l;
velAng += accAng*0.01;
posAng += ((velAng*0.01)-(0.5*accAng*0.0001));
pos[0] = l*sin(posAng);
pos[1] = l*cos(posAng);

These calculations are run every time a timer ticks (which is every 10 milliseconds or 0.01 seconds) and the pos[0] and pos[1] variables are the x and y coordinates of the pendulum bob relative to the origin (the point where the string is held). 'g' is 9.81 and 'l' is 1.5.

For some reason, the pendulum swings noticeably higher on each swing and I can't figure out why. If it helps, all of the variables are defined as type 'double'.

Thankyou for your help!
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
You already increased the velocity before changing the position with the assumption that you still have the old velocity. change the order of the calculations or use variables for old and new state.
 

george7378

DON'T PANIC
Addon Developer
Donator
Joined
Jun 26, 2009
Messages
1,045
Reaction score
0
Points
36
I think the equation for posAng already takes that into account by having a - instead of a + in the middle (it's equivalent to the constant acceleration equation s=vt-0.5at^2 rather than s=ut+0.5at^2).

I'm beginning to think it may be because I'm assuming that the acceleration remains constant for each 10ms timestep when it's actually varying? I'm not exactly sure!
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
I think the equation for posAng already takes that into account by having a - instead of a + in the middle (it's equivalent to the constant acceleration equation s=vt-0.5at^2 rather than s=ut+0.5at^2).

I'm beginning to think it may be because I'm assuming that the acceleration remains constant for each 10ms timestep when it's actually varying? I'm not exactly sure!

You are not understanding the math behind your code there, so in more words:

You start at t0, the beginning of your timestep.

Your state at t0 is v0 and p0.

And the end of the timestep, you want to have t1 (the beginning of the next step) , v1 and p1.

You calculate:
v1 = v0 + a * Δt
p1 = p0 + v1 * Δt - 0.5 * a * Δt²

Instead of:
v1 = v0 + a * Δt
p1 = p0 + v0 * Δt + 0.5 * a * Δt²

See the problem?

Yes, you should in fact also have "a" selected midway between "t0" and "t1" - that is where Runge-Kutta becomes shiny for further improving accuracy.

But currently, alone the fact that you use the wrong sign and the wrong velocity value is enough to explain the increasing energy of your pendulum.
 

george7378

DON'T PANIC
Addon Developer
Donator
Joined
Jun 26, 2009
Messages
1,045
Reaction score
0
Points
36
OK, so I changed it to this:

Code:
[B]double oldVel = velAng;[/B]
accAng = -g*sin(posAng)/l;
velAng += accAng*0.01;
posAng += (([B]oldVel[/B]*0.01)[B]+[/B](0.5*accAng*0.0001));
pos[0] = l*sin(posAng);
pos[1] = l*cos(posAng);

It still does exactly the same thing. I always thought that:

v1 * Δt - 0.5 * a * Δt²
and
v0 * Δt + 0.5 * a * Δt²

were equivalent? Maybe I'm not seeing the problem!

0afdfec75f76bf22c8071ba92f406eca.png


I'm using the 5th equation in this list and you're using the 2nd one.

Maybe if I just stare at the code a bit more it will come to me!
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
The 5th equation is based on the 4th equation.

I will give it a quick try in python.
 

RisingFury

OBSP developer
Addon Developer
Joined
Aug 15, 2008
Messages
6,427
Reaction score
491
Points
173
Location
Among bits and Bytes...
I think this might be purely numerical.

On the uphill, the pendulum calculates acceleration before it increases its altitude. The acceleration there is lower than higher up, so it's less than average. The pendulum gains energy because of that and is able to reach a point much higher up than it would.

The apex is then higher up and the acceleration at that point is also higher. On the down swing, it benefits from the extra bit of energy it accumulated on the uphill. It loses some of it going back down, but I don't think it loses as much as it got getting up.


Try decreasing the time step by 10 or 100 times and run the simulation again. See if that improves the results.
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
Yes, it must be numerical, but shorter time steps would not stop the problem, only reduce its magnitude.

Get the same output in a quick and dirty python script.

I now try it with a RK4 solver, should be better since it also includes acceleration.
 

george7378

DON'T PANIC
Addon Developer
Donator
Joined
Jun 26, 2009
Messages
1,045
Reaction score
0
Points
36
I already tried decreasing the timestep and yes, it did make the problem less noticeable, but it was still there. At the moment it's running at 100 updates per second (the timer ticks every 10ms) and it's still visibly wrong.

I also tried embedding more smaller steps within each 10ms timestep like this:

Code:
for (int i=0; i<10; i++)
{
accAng = -g*sin(posAng)/l;
velAng += accAng*0.001;
posAng += ((velAng*0.001)-(0.5*accAng*0.00001));
}
pos[0] = l*sin(posAng);
pos[1] = l*cos(posAng);

This didn't seem to make it any smoother. I guess the only way round it is to calculate for the middle of each timestep rather than the beginning, like you said...
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
Using RK4 for solving the position and velocity, I get pretty stable calculations now, but I need either a shorter time step or a longer pendulum, so there are less sign-changes along the calculation.

The example looks like this now:
Code:
from math import *
import csv

g = 9.81
l = 1

def a_func(p):
	return -sin(p) * g/l
	
def v_func(p0, v0, dt):
	a0 = a_func(p0)
	
	pa = p0 + v0 * 0.5 * dt
	va = v0 + a0 * 0.5 *dt
	
	pb = pa + 0.5 * dt * va
	vb = va + a_func(pa) * 0.5 * dt
	
	pc = p0 + dt * vb
	vc = v0 + dt * a_func(pc)
	
	
	return (v0 + 2 * va + 2 * vb + vc)/6.0
	
def p_func(p0, v0, dt):
	a0 = a_func(p0)
	
	pa = p0 + v0 * 0.5 * dt
	va = v0 + a0 * 0.5 *dt
	
	pb = pa + 0.5 * dt * va
	vb = va + a_func(pa) * 0.5 * dt
	
	pc = p0 + dt * vb
	vc = v0 + dt * a_func(pc)
	
	
	return (p0 + 2 * pa + 2 * pb + pc)/6.0	
	

	

def pend_func(state, deltaT):
	nstate = state.copy()	
	nstate['t']	+= deltaT	 
	nstate['v'] = v_func(state['p'], state['v'], deltaT);	
	nstate['vdeg'] = degrees(nstate['v'])
	nstate['p'] = fmod(p_func(state['p'], state['v'], deltaT), pi)
	nstate['pdeg'] = degrees(nstate['p'])
	nstate['a'] = a_func(nstate['p'])
	nstate['adeg'] = degrees(nstate['a'])
	nstate['e'] = nstate['p'] * g + 0.5 * nstate['v']**2
	return nstate
	
	
state = {'v':0.0, 't':0.0, 'p':0.25*pi, 'vdeg':0.0, 'pdeg':degrees(0.25*pi), 'a':a_func(0.25*pi), 'adeg':degrees(a_func(0.25*pi))}

csvfile = open('pendulum.csv', 'w', encoding='utf8')
csvw = csv.writer(csvfile, dialect='excel-tab')
csvw.writerow(['t[s]', 'p[°]', 'v[°/s]', 'a[°/s²]'])

while state['t'] < 10.0:
	csvw.writerow([state['t'], state['pdeg'], state['vdeg'], state['adeg']]);
	state = pend_func(state, 0.01)
	
csvfile.close()

After some bugfixing and exporting the samples to csv format, I get a nice sine wave for every parameter, as you can expect, without increases in energy.
 
Last edited:

csanders

Addon Developer
Addon Developer
Joined
Jan 18, 2012
Messages
219
Reaction score
0
Points
0
Location
Plymouth
Nevermind... tab separated...
 
Last edited:

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
Nevermind... tab separated...

Exactly - also the csv.writer handles the delimiter itself, depending on the "dialect" parameter that you specify, no need for a manual delimiter as parameter between the values. I selected Excel simply out of laziness, we use it like that at work often and it reduces the annoyance for synchronizing an Excel sheet with the CSV file data automatically after every bugfix.

I could also have used a DictWriter there.

CSV is "comma separated values" but in cold dark reality it means "a text file with data separated by any characters that is not data" - you must be specify what kind of CSV you actually mean. There is also ";" as delimiter pretty often.
 

RisingFury

OBSP developer
Addon Developer
Joined
Aug 15, 2008
Messages
6,427
Reaction score
491
Points
173
Location
Among bits and Bytes...
I already tried decreasing the timestep and yes, it did make the problem less noticeable, but it was still there.

That's what I wanted to know. That's usually the smoking gun proof that the problem is numerical and not a flaw in your algorithm. If there's a flaw in the algorithm, it would produce the same error (or an error of a higher magnitude than expected) at lower time steps.



I also tried embedding more smaller steps within each 10ms timestep like this:

Code:
for (int i=0; i<10; i++)
{
accAng = -g*sin(posAng)/l;
velAng += accAng*0.001;
posAng += ((velAng*0.001)-(0.5*accAng*0.00001));
}
pos[0] = l*sin(posAng);
pos[1] = l*cos(posAng);

This didn't seem to make it any smoother. I guess the only way round it is to calculate for the middle of each timestep rather than the beginning, like you said...

This produced *exactly* the same result as if you just decreased the time step by 10 times.

Try this:
First calculate where the object would be the same way you're doing now. Then figure out the acceleration at that point. It won't be exactly correct, but close enough. Then average it with the one you got at the beginning of the step and use that average to propagate.

After that, try a Runge-Kutta order 4 as Urvumpe suggested and compare the results of all integrators.
 

csanders

Addon Developer
Addon Developer
Joined
Jan 18, 2012
Messages
219
Reaction score
0
Points
0
Location
Plymouth
I ran the RK4 solution various ways, here were my results for p:

pendulum.png
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
37,588
Reaction score
2,312
Points
203
Location
Wolfsburg
Preferred Pronouns
Sire
:rofl::lol: Did you put a wall 5 seconds into the simulation? It happens at roughly t = 5 s on every graph. Why?

No, he just post-processed the 100s data to show the first and last 5 seconds at different sampling rate. :thumbup:
 

csanders

Addon Developer
Addon Developer
Joined
Jan 18, 2012
Messages
219
Reaction score
0
Points
0
Location
Plymouth
:rofl::lol: Did you put a wall 5 seconds into the simulation? It happens at roughly t = 5 s on every graph. Why?

1) printed the first 5 seconds and last 5 seconds to file, to better see the accumulation error.
2) opened file in excel and made graph from printed data

PHP:
while state['t'] < 100.0:
	if state['t'] < 5 or state['t'] > 95:
		csvfile.writerow(state['t'], state['pdeg'], state['vdeg'], state['adeg'])
	state = pend_func(state, 0.1)
 
Last edited:

george7378

DON'T PANIC
Addon Developer
Donator
Joined
Jun 26, 2009
Messages
1,045
Reaction score
0
Points
36
Thanks for the pointers everyone - I did the average of the acceleration at the beginning and the end of each timestep to do the propagation and it looks like it's working now! In fact, here's a video of how it looks:


I also did a double pendulum as you can see, but I didn't use the averaging in this one - it's far more chaotic so it's harder to notice the small inaccuracies in momentum.
 

george7378

DON'T PANIC
Addon Developer
Donator
Joined
Jun 26, 2009
Messages
1,045
Reaction score
0
Points
36
Thanks :) It always amazes me to see equations come to life like that!
 
Top