数学研发论坛

 找回密码
 欢迎注册
楼主: KeyTo9_Fans

[原创] 谢尔宾斯基地毯与圆的交集

  [复制链接]
发表于 2017-9-25 14:09:56 | 显示全部楼层
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$

顺便贴一下我画的图


点评

现在好了,7次迭代后,8^7个格子,圆内的有1543172个,占比0.735842  发表于 2017-9-25 22:00
结果不正确吧,应该在3/4附近  发表于 2017-9-25 15:49
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-25 16:53:01 | 显示全部楼层
{3,17/64},{4,33/128},
显然不对,圆内的应该只增不减
另外最好再给出圆上的数目,这样就可以非常容易的同时算出上下界
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-25 17:32:51 | 显示全部楼层
mathe 发表于 2017-9-25 16:53
……………………
最好再给出圆上的数目,这样就可以非常容易的同时算出上下界


取上下界的中值,即圆周上的非白方块只计一半,收敛速度会快得多吧,很早就能看到大致结果。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-25 21:43:17 | 显示全部楼层
老大们说的对,重新计算了下:
  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}
复制代码

以不同的颜色标记,画图如下:

Untitled-1.png

Untitled-2.png

点评

Norm完Sort一下,取第一个和最后一个  发表于 2017-9-26 08:39
给定一个格子的中心点,求出格子的四个顶点中,离圆心{0,0}最近和最远的顶点的坐标  发表于 2017-9-26 08:32
就是针对一个序列,在一趟计算中 同时求出其中的最大值和最小值  发表于 2017-9-26 08:31
MinMax[], 俺M10中没有这个函数,是否可用Min@Max代替?  发表于 2017-9-26 08:28
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-26 00:46:41 | 显示全部楼层
@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 你是按小方块的中心来划分的么?

点评

软件难免有bug的,一般还不止一处呢,慢慢来,三番五次就消除了。  发表于 2017-9-26 08:33
根据中心点获取四周的顶点的时候乘以上一个边长的1/6,或者当前边长的1/2,我搞错位了。代码已经更正  发表于 2017-9-26 08:15
惭愧,只需要将其中一处/2改成/6就行了,已经更正  发表于 2017-9-26 07:52
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-26 06:48:53 来自手机 | 显示全部楼层
看wayne给出的图,蓝色圆(边界上的)太多了

点评

根据中心点获取四周的顶点的时候乘以上一个边长的1/6,或者当前边长的1/2,我搞错了。代码已经更正  发表于 2017-9-26 07:55
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-26 06:56:25 来自手机 | 显示全部楼层
这道题目还可以用Bresenham算法求出和圆相交各小正方形。而对于每个这样小正方形,实际上可以使用它两坐标的三进制形式计算出其下圆内黑格数目。这样做的好处是空间复杂度比较小
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-26 14:25:55 | 显示全部楼层
这个题目还可以用蒙特卡洛法,然后我计算划分到20层后的蒙特卡罗法,结果在0.736左右

点评

蒙特卡洛法精度有限,我这个是跑1000000次,精度大概在1/1000左右,所以这个结果属于正常范围  发表于 2017-9-29 09:14
第20层的结果是 212923784529434167/288230376151711744 = 0.73872777523407348252  发表于 2017-9-28 14:59
大于0.73871530  发表于 2017-9-26 22:19
根据单调性,好像大于 0.73769951  发表于 2017-9-26 20:44
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-26 18:26:38 | 显示全部楼层
本帖最后由 lsr314 于 2017-9-26 18:29 编辑

不考虑与圆相交而被切割的误差,只以小正方形的中心位置来划分圆内圆外,那么可以用三进制来判断,代码如下:
  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}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-26 20:26:35 | 显示全部楼层
本帖最后由 chyanog 于 2017-9-26 20:35 编辑

写了两个版本的Mathematica代码:
第一个需要11.1版本,使用了MengerMesh和计算几何相关的函数,代码简单易懂,不过效率不高
  1. Clear["`*"];

  2. n=3;

  3. poly=MeshPrimitives[MengerMesh[n],2];

  4. Graphics[{
  5. Antialiasing->False,
  6. EdgeForm[Black],
  7. Red,Select[poly,RegionWithin[Disk[{0.5,0.5},0.5],#]&],
  8. Blue,Select[poly,!RegionDisjoint[Circle[{0.5,0.5},0.5],#]&],
  9. Gray,Select[poly,!RegionWithin[Disk[{0.5,0.5},0.5],#]&&RegionDisjoint[Circle[{0.5,0.5},0.5],#]&],
  10. Thick,Magenta,Circle[{0.5,0.5},0.5]
  11. }]
复制代码

第二个版本速度快一些,不限于版本11
  1. Clear["`*"];

  2. carpet[n_]:=Nest[ArrayFlatten[{{#,#,#},{#,0,#},{#,#,#}}]&,1,n];

  3. n=3;

  4. ArrayPlot[carpet[n],Mesh->All]

  5. p=Rescale[{{#1,#2},{#1+1,#2},{#1+1,#2+1},{#1,#2+1}}&@@@N@Position[carpet[n],1]];

  6. Graphics[{
  7. EdgeForm[Black],
  8. Antialiasing->False,
  9. Red,Polygon@Select[p,AllTrue[#,Norm[#-0.5]<0.5&]&],
  10. Blue,Polygon@Select[p,With[{m=Norm[#-0.5]&/@#},Min@m<0.5&&Max@m>0.5]&],
  11. Gray,Polygon@Select[p,AllTrue[#,Norm[#-0.5]>0.5&]&],
  12. Thick,Magenta,Circle[{0.5,0.5},0.5]
  13. }
  14. ]
复制代码




补充内容 (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}
2017-09-26_202847.png

点评

光画图,不计算一下“面积”么?  发表于 2017-9-26 21:02
赞  发表于 2017-9-26 21:02
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-12-10 19:24 , Processed in 0.071367 second(s), 18 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表