I'll try to give you a background on where the problem occured. The problem itself is defined near the end of the post.

I'm running a sort of a simulation of a particle in a box. I finally got my hands on Mathematica10 so I decided to use new built-in mechanisms such as regions. Long story short: I need to determine the point at which the particle hits the bounding region.

I defined my box by:

{X, Y, Z} = {.05, .05, .3};
ip1 = InfinitePlane[{{0, 0, 0}, {0, Y, 0}, {0, Y, Z}}];
ip2 = InfinitePlane[{{X, 0, 0}, {X, Y, 0}, {X, Y, Z}}];
ip3 = InfinitePlane[{{0, 0, 0}, {X, 0, 0}, {X, 0, Z}}];
ip4 = InfinitePlane[{{0, Y, 0}, {X, Y, 0}, {X, Y, Z}}];
ip5 = InfinitePlane[{{0, 0, 0}, {0, Y, 0}, {X, Y, 0}}];
ip6 = InfinitePlane[{{0, 0, Z}, {0, Y, Z}, {X, Y, Z}}];
reg = RegionIntersection[Cuboid[{0, 0, 0}, {X, Y, Z}], RegionUnion[ip1, ip2, ip3, ip4, ip5, ip6]];

I generate the initial position and velocity of the particle as random numbers, then draw a HalfLine from the position of particle along the velocity vector. Finally, I check the intersection of the HalfLine and the bounding region reg with RegionIntersection.

Then I can extract the coordinates of the bouncing point and calculate the new velocity vector, define a new HalfLine, find the intersection, and so on...

The Problem

At some point, however, I get an error caused by the point not being a member (or element, both functions do the same in this case) of my bounding region (sounds ridiculously, because it's supposed to be the intersection of the region and the halfline).

When I print the coordinates of that point, I get:

In: p1
Out: {0.0000316867, 0.00969779, 0.3}

So everything looks fine, but doesn't work. I checked it with Manipulate - I entered these coordinates manually and RegionMember[reg,{0.0000316867, 0.00969779, 0.3}] evaluated to True.

Then I copied the coordinates of p1 from the output above and found out that they were in fact:

{0.00003168673646693154`, 0.009697790214128308`, 0.29999999999999993`}

The last coordinate was at fault. Is this a bug in the new functions? Any ideas how to get around it?

EDIT: Computing p1: @user21

line = HalfLine[r + v/10^8, v];
p1 = RegionBounds[RegionIntersection[line, reg]][[All, 1]];

If line would start on the wall, RegionIntersection would consist of two points. To avoid that, I translate the origin of the HalfLine along v a bit.

EDIT 2: I've found a way around, but it doesn't answer the question about precision

Okay, instead of defining p1 as an element of RegionBounds I used RegionCentroid:

old: p1 = RegionBounds[RegionIntersection[line, reg]][[All, 1]];
new: p1 = RegionCentroid[RegionIntersection[line, reg]]

And it seems to have solved the problem. But I have no idea why, because for a single point both RegionCentroid and RegionBounds should give the same result.

I tried both functions on a "buggy" (i.e. found by RegionBounds) point and here's what I got:

In: RegionBounds[RegionIntersection[line, reg]]
Out: {{0.0494424, 0.0494424}, {0.0341449, 0.0341449}, {0.3, 0.3}}

And after copy-pasting the output:

{{0.04944240524519992`, 0.04944240524519992`}, {0.03414494309476567`, 0.03414494309476567`}, {0.30000000000000004`, 0.30000000000000004`}}

And for RegionCentroid:

In: RegionCentroid[RegionIntersection[line, reg]]
Out: {0.0494424, 0.0341449, 0.3}

Copy & paste:

{0.04944240524519991`, 0.03414494309476552`, 0.3`}

Well, I can imagine that it's due to internal functions used, but I'm still somehow confused that the results do not "converge" for this case.

  • $\begingroup$ how did you compute p1? $\endgroup$ – user21 Aug 1 '14 at 14:34
  • $\begingroup$ @user21 I added this information to the post. $\endgroup$ – Daroslav Aug 1 '14 at 14:41
  • $\begingroup$ What are r and v? You could use exact arithmetic with Rationalize[{.05, .05, .3}] and making sure r and v are exact too. $\endgroup$ – user21 Aug 4 '14 at 12:15
  • $\begingroup$ r and v are the coords and the velocity of my particle, respectively. I will try that and see what happens. Thank you! $\endgroup$ – Daroslav Aug 4 '14 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.