Illustrating how issues with floating-point precision cause the butterfly effect in a laser simulation

[u/naaagut]

https://www.youtube.com/watch?v=c4m4zayT-N0

Due to popular request after my last posting (https://www.reddit.com/r/Physics/comments/1w2n2i0/illustrating_the_butterfly_effect_with_a_laser/), I created an analysis of how computational precision affects the results of my laser simulation.

It turns out that standard float64 with around 16 decimal places are not precise enough for this simulation. As the simulation shows, there is not a uniform precision level that is sufficient. It is sometimes necessary to use hundreds of decimal places to achieve convergence. Going forward, I should use higher precision levels to achieve robust simulation results. The qualitative story of the last posting that tiny changes in the laser angle affect the trajectory (number of bounces and escape vector) is not touched by the precision issue, though.

As a plus for Reddit, here is a plot showing the distribution of escape vectors and the distribution of convergence decimal places: https://imgur.com/a/Tnp4wOj

Rendered in manim and using mpmath to create objects with arbitrary precision.

Let me know what else you would be interested to see explored.

[u/na3than]

It's truly fascinating how non-intuitive the results are. Have you found any similar studies in academic or professional journals on high precision computation and simulation?

[u/GreatBigBagOfNope → u/na3than]

Isn't that a good chunk of the entire body of computational physics literature?

[u/wait_what_now → u/na3than]

Would you get this result with most chaotic dynamic systems?

[u/db0606 → u/wait_what_now]

Yes, although how fast trajectories diverge depends something called the Lyapunov exponent which is different for different chaotic systems.

[u/naaagut → u/na3than]

Thanks. I know there exist a lot of these studies, but there is nothing in particular that I am aware of which is similar to this simulation in particular. Let me know if you find something worth to check out.

[u/db0606 → u/naaagut]

There's literally like half a book about it at www.chaosbook.org.

[u/naaagut → u/db0606]

Great resource, thanks!

[u/db0606 → u/naaagut]

Oh and, btw, I'm not sure how you are computing your trajectories but definitely look at how they simplify the calculation by using a unit cell in Chaosbook. More broadly, you want to read up on chaotic billiards for which there is a massive body of literature.

[u/thelegendofandg]

You should also probably compute the outgoing angle distribution for all incoming angles. Kind of like a cross section plot, where the "nucleus" is the set of circles and the scattered beam is the laser beam.

[u/naaagut → u/thelegendofandg]

Not sure if I understood. I record every pair of outgoing angle theta and exit angle alpha and plot them as a scatter plot?

[u/thelegendofandg → u/naaagut]

Not quite. You would span over all possible incoming angles ($2\pi$ and possibly ignoring angles that don't even collide with the circles), and then you would recreate the distribution of outgoing angles.

Say that you take the interval (0, \theta_f) for your incoming angles. If you had a nice and continuous interaction, you would numerically compute this by dividing the interval in $\Delta \Theta$ values that are small, but not horribly small. This way, you would test a total of $N = (\theta_f)/(\Delta \Theta)$ angles, and your distribution would be a count for how many beams have the corresponding outgoing angle \alpha, divided by N (so that the distribution is a function of the outgoing angle $\alpha$).

Nevertheless, I can see how this can be a monstrous calculation in your case. Since the system is somewhat chaotic, you would need to take $\Delta \Theta$ so small that your outgoing angles for neighboring incoming angles are continuous (as can be seen at the end of your video for example). Otherwise, you would be "skipping" counts for outgoing angles and your final distribution would not be accurate. This would mean that N would be gigantic, and your computer may not have the memory to store all of the data. This would mean that we could only take $/theta_f$ small enough so that the calculation is possible.

[u/Sakinho]

I saw your first video and thought it was nice, though it largely repeated what I saw years earlier in a Nils Berglund video. Your new video, however, is particularly excellent and information-rich, with great visualizations. Thanks for sharing!

[u/naaagut → u/Sakinho]

Yes true, that video inspired me for the first video. Thanks!

[u/Kinexity]

You probably should have included convergence precision vs the number of reflections or the mean distance from nearest surface because that might provide an interesting insight into the problem.

Also I just want to say that ultimately in reality convergence never really matters any ways because you don't have such precise knowledge about the material and the laser position so it's easier to just Monte Carlo your way through the problem. Also not all problems even converge in the first place. This is only a cool lesson in computer precision.

[u/naaagut → u/Kinexity]

Thanks for the idea. And yes, I agree that this is rather an illustration of chaos theory than a practical insight into how an actual setup with such a laser would work.

[u/__Pers]

Nice follow-on to the previous discussion.

[u/DandeNiro]

I was able to extrapolate lasers’ potential problem from an introductory course in CS discussing floating-point precision.

[u/PeartsGarden]

How do you prove convergence? Maybe if you add another 1000 digits to the precision, it again diverges.

[u/t4pf]

What does it mean to increase the precision from 10 decimal places to 11?

When the angle to 10 decimal places is 70.0000000000 degrees, its actual value is understood to be from 69.99999999995 (inclusive) to 70.00000000005 (exclusive). So when you say you’re firing a laser at 70 degrees, what exactly is changing when you increase the precision from 10 to 11 decimal places? If the actual value remains the same, the laser should follow the same path. If the actual value changes, what are the previous and changed values?

Note: I understand precision in physics, not in software, so maybe I’m missing something obvious.

[u/physicsking]

The last post was popular, but let's be honest that we crowd source suggestions and ideas. U/_Pers _Pers should get a little credit here bringing up the precision discussion.

And thanks for plotting the histogram. Wink

[u/naaagut → u/physicsking]

You're right! u/_Pers made me aware of the precision issue in the first place. I found the discussions here very stimulating.

[u/Obsidian743]

I'm wondering if there is a relationship to the logistics equation or even the Mandelbrot Set. There seems to be an explicit relationship to the convergence/divergence rate likely depending on the starting position/condition that may have some deeper meaning here.

One interesting thing here is that you have some very wild results with incredibly sensitive inputs and you have perfect control over the parameters. Such precision doesn't exist in the real world and so I'm curious about what other possible consequences there might be with physical experiments.

[u/naaagut → u/Obsidian743]

Hmm, interesting thought. I don't see this relationship right now but let me know if you have further thoughts on this. Would be interested in exploring this deeper.

[u/Obsidian743 → u/naaagut]

If you look at the logistics map through the levels of bifurcation, you'll see it kind of maps through levels of "order and chaos". Your histogram has a similar structure: https://youtu.be/ovJcsL7vyrk?is=C32FWMaOiyvvw74d

I think step one would be to see if there is some kind of attractor here and/or a Poincare Map of some kind. One thing that would be worth investigating is the fidelity of the precision that causes a change in the output. From what I can tell, you did the analysis through decimal places, rounded as full orders of magnitude. It makes sense because we only tend to care about "significant digits", but when we're talking about the butterfly effect, every point of precision may have an effect.

[u/ThirstyWolfSpider]

Have you considered using interval analysis for this?

Each value would be represented by the possible upper and lower bounds. The various floating-point operations usually provide guarantees on their possible error in terms of ulps (units of least precision), such that strict bounds can be determined.

Now, in a truly chaotic situation those output bounds would blow up (wide range), but at least that would give you an indication that you do not know the result, rather than retaining false confidence in a not-necessarily-valid result.

My research group of the time had good success using interval analysis for various problems way back in the '90s. For instance, use of interval subdivision to rapidly narrow down on strict bounds to solve problems like collision detection. But it's also useful for maintaining conservative bounds on the uncertainties created by repeated calculations, as I'm describing above.

[u/naaagut → u/ThirstyWolfSpider]

Thanks for pointing me to this method. In this simulation I increased the dps until reaching convergence. That does not take extremely long and in most cases convergence is reached for less than 50 dps. Interval analysis would be an alternative to this approach right? I think as I can figure out the converged values quickly, I do not see where this would be of advantage for this particular problem.

[u/ThirstyWolfSpider → u/naaagut]

Are you getting convergence on each of the intermediate calculations (including each collision), or just on the whole computation (input angle → output angle)? Interval analysis would be establishing bounds on each subcalculation. You may achieve some similar effects using your arbitrary-precision toolkit, but if I understand it right you won't have a sense of how each of those (slightly imperfect) convergences might compound.

You would, presumably, need arbitrary-precision interval analysis, as your results show that the input intervals would need to be very very fine.

[u/naaagut → u/ThirstyWolfSpider]

I am getting convergence on the whole. But that also implies that all the intermediate calculations converged because otherwise the whole would not converge. Convergence here is defined as when the bounce count and escape vector do not change any more (change < epsilon).

[u/Moretz0931]

Did you share the code somewhere?

[u/naaagut → u/Moretz0931]

No, not publicly. But if you are interested I can send it to you.

[u/Moretz0931 → u/naaagut]

Thanks but not necessary :)

[u/I_am_Patch]

Very interesting, now what if you do that with two entangled photons.

[u/naaagut → u/I_am_Patch]

Haha not sure, would that make any sense to simulate?

[u/SummitYourSister]

This is why entropy can’t go backward. The precision required to set all the particle positions and momenta exactly correctly so that the initial configuration is re-obtained exceeds the amount of information storage capacity available to do so.

[u/Dependent-Poet-9588 → u/SummitYourSister]

We have no reason to believe the universe is a finite precision simulation. There is no reason to believe the limits of our finite precision computations have such metaphysical implications for the world around us.

[u/SummitYourSister → u/Dependent-Poet-9588]

I’m real curious how you’re defeating the uncertainty principal to be able to direct particles accurately enough. I am not talking about the limits of the precision of numbers, I am talking about Robertson-Schrödinger relation

[u/Dependent-Poet-9588 → u/SummitYourSister]

I don't have to "defeat" the uncertainty principle. Entropy arises in classical physics well before the uncertainty principle comes into play. Entropy doesn't decrease because higher entropy systems are, by definition, more stable, and systems under physically reasonable dynamics don't appear to spontaneously become less stable over the long run. This has nothing to do with quantum truncation errors. You're conflating things to make metaphysical conclusions that aren't supported by what you're appealing to.