How can I estimate pi by using a unit hexagon in Matlab?

We can do easily on estimate pi using 1-by-1 square tiling, where codes are below. But how about “hexagons”? Still use for loop?

% Pi Via Tiling

% Enter the disk radius...
clc
n = input('Enter an  integral radius n: ');

% Tile just the first quadrant, then multiply by four...
N1 = 0;
for k = 1:n
    % Add in the number of uncut tiles in row k...
    m = floor(sqrt(n^2 - k^2));
    N1 = N1 + m;
end

% Display the estimate...
rho_n = 4*N1/n^2;
clc
fprintf('n      = %d\n',n)
fprintf('rho_n  = %12.8f\n',rho_n)
fprintf('Error  = %12.8f\n',abs(pi-rho_n))

[Image Analyst]

The problem statement doesn't mention pi. So why do you want to measure it?

[Zian Fu]

Yes, I wanna estimate pi using hexagons tiling method, maybe the same thinking as "squares", but I have no idea how to count the number of hexagons which I guess, is crucial.

Answer 1

[John D'Errico]

I won't do your homework for you. So you need to think. I'll give you a few (huge) hints:

First, is there some symmetry you can take advantage of? It won't kill you if you cannot see this, but symmetries are often valuable tools to employ. They reduce the programming effort, and the time to solve a problem.

Can you define the coordinates of the center of any tile? The homework even gave you a big hint there, in the form of * versus o tiles. Come on. This is basic geometry.

If you know the center of any tile, then can you find the coordinates of the 6 corners of the hex? Again, basic geometry, plus a little trig. Start with a hex centered at (0,0). Are all the hexes of fundamentally the same shape and size? Is there some simple, additive operation you can do?

If a hex is FULLY inside the circle then you will count it. How can you tell if a hex is inside a circle? Is there a simple scheme you can find, based on the vertices of the hex in question?

When you see a problem too big and difficult for you to get your head around, break it into small, bite sized chunks. Make big problems into small ones, then solve the small problems, one at a time. Put it all together, and you are done. This advice applies to any problem you will be forced to solve in the future.

Or, if all else fails, you can just do this:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199

Actually, that last digit was rounded off. It would take REALLY small hexagons to get that degree of accuracy though. :)

[Zian Fu]

Maybe it's a little hard for me to find a suitable coordinate with * and o. Can you give me more hints?

[John D'Errico]

Sorry. I don't like to provide code for homework problems, and if I gave you much more, then I'd be writing your homework.

THINK. Basic geometry. How big is a hexagon, from side to parallel side? Assume that you know the side length of one hexagon. That tells you everything about the tiling.

If you know the width of a hex from side to side, then what is the distance between the centers of two adjacent hexagons that share a face? This applies for each of the * and o hexagons. A pair of adjacent hexagons can be viewed for example, as sharing the same x coordinate, with different y.

[Image Analyst]

An array of hexagon centers is simply two superimposed rectangular grids, one halfway to the right and upwards of the other, as I showed in my answer.

[Zian Fu]

Firstly, it's not my homework.:)

Secondly, here is my code below, which, I guess, is too complex because of too many "for" increasing the time and somewhere isn't very thoughtful, for example, the expression for ending value needn't reach "r" (a number less than r may be enough), though I think it won't affect the correct final result.

I found an unacceptable error after running my program, that difference between pi and the estimated is 0.08 even if radius=100, which is much worse than squares' estimation. And then I improved my code somewhere, but the significant improvement still didn't appear.

clc;
clear;
tic
r=input('Enter a real numberas radius=');
count_star=0;
count_circle=0;
% Due to symmetry, consider 1/4, upper right area
% Consider "*" tiles firstly
tic
for x_star= 0:3:r %loop for coordinate x of "*" tile center
    for y_star= sqrt(3)/2:sqrt(3):r %locate heaxgon's upper right corner
        if x_star+1<=r && (x_star+0.5)^2+y_star^2<=r^2
            %guarantee hexagon's upperright and rightmost is inside circle
            count_star=count_star+1;
            %symmetry can lead to more heagons counted along positive x 
            %and y axis, so we have to extract them specially
            if x_star==0
                eliminate_ystar=count_star;
                %in this way, we can get the final numbers along positive
                %y axis by loop
            end
            if y_star==sqrt(3)/2
                eliminate_xstar=x_star;
                %along positive x asix by loop
            end
        end
    end
end
toc

tic
% Consider "o" tiles secondly
for x_circle=2:3:r %condinate x of "o" tile's lower right corner
    for y_circle=sqrt(3):sqrt(3):r % Similarly, y of upper right corner
        if x_circle+0.5<=r && x_circle^2+y_circle^2<=r^2
            % guarantee hexagon's upperright and rightmost is inside circle
            count_circle=count_circle+1;
        end
    end
end
toc
% Make use of symmetry, calculate the total
count_total=4*(count_circle+count_star);
% eliminate the repetitive cases which can only be made by "*" tiles
% the last "1" is the origin
count_total=count_total-2*(eliminate_ystar-1)-2*(eliminate_xstar-1)+1;
% Display the estimate and relative errors
S_total=(3*sqrt(3)/2)*count_total; % Total area of hexagons
estimate_pi=S_total/(r^2);
fprintf('radius    %30.8f\n',r)
fprintf('estimate_pi   %22.8f\n',estimate_pi)
fprintf('Error      %25.8f\n',abs(pi-estimate_pi))
toc

So what I really want to ask are:

1. What's the biggest problem in my code that leads to the unacceptable error?

2. Except that are there any details that can be improved to save run time or decrease the complexity?

Thank you very much!

[John D'Errico]

Is this homework? :) It still sort of is homework, IF you are using it to learn MATLAB, or anything about mathematics. You gain by the doing. I never have a problem in helping someone who is making an effort though, someone who wants to learn.

I'll come back to this later when I have more time to read through it. Must run out this am though.

The one thing I saw immediately is this:

if y_star==sqrt(3)/2

A BAD idea. NEVER test for exact equality of floating point numbers like that. Use a tolerance instead. Perhaps this:

if abs(y_star - sqrt(3)/2) < 10*eps

[Zian Fu]

Thank you. waiting for your more specific reply:)

Answer 2

[Image Analyst]

A little bit of code adapted from the FAQ to get you started:

r = 11.1234; % Whatever.
X = [0 : 3 : ceil(r), 0 : -3 : floor(-r)];
Y = [0 : 3 : ceil(r), 0 : -3 : floor(-r)];
[x1, y1] = meshgrid(X, Y);
x2 = x1 + 1.5;
y2 = y1 + 1.5;
x = [x1(:), x2(:)];
y = [y1(:), y2(:)];
plot(x, y, 'b.', 'MarkerSize', 20);
grid on;
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';

% Plot the circle.
centerX = 0;
centerY = 0;
hold on;
rectangle('Position',[centerX - r, centerY - r, r*2, r*2],...
    'Curvature',[1,1],...
  'EdgeColor', 'r',...
  'LineWidth', 2,...
    'FaceColor','none');
axis square;

See if you can figure out, from the first code chunk in the FAQ: http://matlab.wikia.com/wiki/FAQ#How_do_I_create_a_circle.3F how to count how many of the (x,y) points are inside of the circle.

Use input or inputdlg() to ask the user for the radius.

[John D'Errico]

A requirement was that a hexagon needs to be ENTIRELY inside the circle to be counted. So be careful.

[Image Analyst]

Yes, it's a bit tricky. I didn't show the hexagons themselves, just their centers, and you can't just see if the centers are within the circle. You have to see if any little tip of the hexagon extends outside of the circle. There will be math ahead, so be prepared.