谢尔宾斯基地毯与圆的交集

wayne

wayne 发表于 2017-9-24 23:59
先占个坑。
对于给定分辨率的“孔方兄”，一个像素格子可以由其中心点的坐标以及边长决定。标记成三元组{a ...

计算出来每迭代一次，小正方形的个数乘以8,最终格子中心落在圆内的个数是  {{1,0},{2,12},{3,136},{4,1056},{5,8536},{6,68488},{7,547928}}
除以对应的总数，得到圆内的格子的占比比例，{{1,0},{2,3/16},{3,17/64},{4,33/128},{5,1067/4096},{6,8561/32768},{7,68491/262144}}
数值结果就是  {{1.,0.},{2.,0.1875},{3.,0.265625},{4.,0.257813},{5.,0.260498},{6.,0.261261},{7.,0.261272}}

极限值看上去像是$ 0.261272$

顺便贴一下我画的图

mathe

{3,17/64},{4,33/128},
显然不对，圆内的应该只增不减
另外最好再给出圆上的数目，这样就可以非常容易的同时算出上下界

hujunhua

mathe 发表于 2017-9-25 16:53
……………………
最好再给出圆上的数目，这样就可以非常容易的同时算出上下界

取上下界的中值，即圆周上的非白方块只计一半，收敛速度会快得多吧，很早就能看到大致结果。

wayne

老大们说的对，重新计算了下：

1. kernel=DeleteCases[Tuples[{-1,0,1},{2}],{0,0}];
   
2. kernel={{-1,-1},{0,-1},{1,-1},{1,0},{1,1},{0,1},{-1,1},{-1,0}};
   
3. rect={{-1,-1},{-1,1},{1,1},{1,-1}};
   
4. g[points_]:=Module[{ps=points},{Flatten[Table[Table[Join[k,{t=MinMax[Norm[ps[[2]]/6 #+k]&/@rect];(Sign[1/2-t[[1]]]+Sign[1/2-t[[2]]])/2}],{k,ps[[2]]/3 #+p[[1;;2]]&/@kernel}],{p,ps[[1]]}],1],ps[[2]]/3}]
   
5. pp=Nest[g,{{{0,0,1}},1},3];
   
6. Graphics[{EdgeForm[Thickness[.001]],Flatten[Table[{Which[j[[3]]==-1,Gray,j[[3]]==0,Blue,j[[3]]==1,Red],Rectangle[j[[1;;2]]+{-1,-1}*pp[[2]]/2,j[[1;;2]]+{1,1}*pp[[2]]/2]},{j,pp[[1]]}],1],Thickness[.005],RGBColor[1,0,1],Circle[{0,0},1/2]}]

复制代码

统计所有的圆外，圆上，圆内的格子数目分别如下：

1. {{1},{1.},8}
   
2. {{1/16,7/16,1/2},{0.0625,0.4375,0.5},64}
   
3. {{13/64,19/128,83/128},{0.203125,0.148438,0.648438},512}
   
4. {{123/512,51/1024,727/1024},{0.240234,0.0498047,0.709961},4096}
   
5. {{259/1024,145/8192,5975/8192},{0.25293,0.0177002,0.72937},32768}
   
6. {{2117/8192,389/65536,48211/65536},{0.258423,0.00593567,0.735641},262144}
   
7. {{34119/131072,1045/524288,386767/524288},{0.26030731,0.0019931793,0.73769951},2097152}
   
8. {{547243/2097152,2801/4194304,3097017/4194304},{0.26094580,0.00066781044,0.73838639},16777216}
   

复制代码

以不同的颜色标记，画图如下：

hujunhua

@wayne 你的结果还是有问题吧。
所谓圆内/外，指的是整个小方块完全在圆内/外，剩下的就是被圆周所穿过的方块。
按｛迭代次数，{圆上，圆内}，{圆内+圆上/2,"面积"}}来陈列，前 3 项数据人工也容易数得出，为
{1, {8, 0}, {4, 1/2}}
{2, {28, 32}, {46, 23/32=0.71875}}
{3, {76, 332}, {370, 185/256=0.72265625}}
第4次迭代
{4, {204,2908}, {3010, 1505/2048=0.73486328125}}

@wayne 你是按小方块的中心来划分的么？

mathe

看wayne给出的图，蓝色圆(边界上的)太多了

mathe

这道题目还可以用Bresenham算法求出和圆相交各小正方形。而对于每个这样小正方形，实际上可以使用它两坐标的三进制形式计算出其下圆内黑格数目。这样做的好处是空间复杂度比较小

mathe

这个题目还可以用蒙特卡洛法，然后我计算划分到20层后的蒙特卡罗法，结果在0.736左右

lsr314

不考虑与圆相交而被切割的误差，只以小正方形的中心位置来划分圆内圆外，那么可以用三进制来判断，代码如下：

1. f[x_, y_] :=
   
2. MemberQ[(IntegerDigits[Min[x, y], 3] - 1)^2 +
   
3.    Take[IntegerDigits[Max[x, y], 3] -
   
4.       1, -Length[IntegerDigits[Min[x, y], 3]]]^2, 0]
   
5. ff[l_] := (s = t = 0;
   
6.   Do[If[f[i, j], ,
   
7.     If[((i + 1/2)/3^l - 1/2)^2 + ((j + 1/2)/3^l - 1/2)^2 < 1/4,
   
8.      s = s + 1, t = t + 1]], {j, 0, 3^l - 1}, {i, 0, 3^l - 1}];
   
9.   s/(s + t))
   
10. Do[Print[{l, N[ff[l]]}], {l, 7}]

复制代码

结果如下：
{1,1.}
{2,0.8125}
{3,0.734375}
{4,0.742188}
{5,0.739502}
{6,0.738739}
{7,0.738728}

chyanog

写了两个版本的Mathematica代码：
第一个需要11.1版本，使用了MengerMesh和计算几何相关的函数，代码简单易懂，不过效率不高

1. Clear["`*"];
   
2. 
   
3. n=3;
   
4. 
   
5. poly=MeshPrimitives[MengerMesh[n],2];
   
6. 
   
7. Graphics[{
   
8. Antialiasing->False,
   
9. EdgeForm[Black],
   
10. Red,Select[poly,RegionWithin[Disk[{0.5,0.5},0.5],#]&],
    
11. Blue,Select[poly,!RegionDisjoint[Circle[{0.5,0.5},0.5],#]&],
    
12. Gray,Select[poly,!RegionWithin[Disk[{0.5,0.5},0.5],#]&&RegionDisjoint[Circle[{0.5,0.5},0.5],#]&],
    
13. Thick,Magenta,Circle[{0.5,0.5},0.5]
    
14. }]

复制代码

第二个版本速度快一些，不限于版本11

1. Clear["`*"];
   
2. 
   
3. carpet[n_]:=Nest[ArrayFlatten[{{#,#,#},{#,0,#},{#,#,#}}]&,1,n];
   
4. 
   
5. n=3;
   
6. 
   
7. ArrayPlot[carpet[n],Mesh->All]
   
8. 
   
9. p=Rescale[{{#1,#2},{#1+1,#2},{#1+1,#2+1},{#1,#2+1}}&@@@N@Position[carpet[n],1]];
   
10. 
    
11. Graphics[{
    
12. EdgeForm[Black],
    
13. Antialiasing->False,
    
14. Red,Polygon@Select[p,AllTrue[#,Norm[#-0.5]<0.5&]&],
    
15. Blue,Polygon@Select[p,With[{m=Norm[#-0.5]&/@#},Min@m<0.5&&Max@m>0.5]&],
    
16. Gray,Polygon@Select[p,AllTrue[#,Norm[#-0.5]>0.5&]&],
    
17. Thick,Magenta,Circle[{0.5,0.5},0.5]
    
18. }
    
19. ]

复制代码

补充内容 (2017-9-27 10:46):
次数，圆上，圆内，圆外，“面积”
{1, 8, 0, 0, 0.5}
{2, 28, 32, 4, 0.71875}
{3, 76, 332, 104, 0.72265625}
{4, 204, 2908, 984, 0.73486328125}
{5, 580, 23900, 8288, 0.73822021484375}