免費論壇 繁體 | 簡體
Sclub交友聊天~加入聊天室當版主
分享
返回列表 发帖

以正奇数边形的两顶点作的等腰三角形的顶点的共线性

本帖最后由 hbghlyj 于 2021-1-14 23:44 编辑
  1. n=27;
复制代码
先把不相邻的顶点的配对存储到pair,它有n(n-3)/2个元素.
  1. points = Table[{Cos[2 k Pi/n], Sin[2 k Pi/n], k}, {k, 0,
  2.    n - 1}]; pair = Subsets[points, {2}];
复制代码
定义以(x1,y1),(x2,y2)为圆心,r为半径的两圆的交点的函数f,当没有交点时输出Nothing(这个对象放在列表中会自动被删除).
  1. f[{{x1_,y1_,k1_},{x2_,y2_,k2_}}]=If[(x1-x2)^2+(y1-y2)^2<=4r^2,{{((x1+x2)((x1-x2)^2+(y1-y2)^2)-Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)])/(2((x1-x2)^2+(y1-y2)^2)),1/(2((x1-x2)^2+(y1-y2)^2)(y1-y2))((x1-x2)^2 y1^2+y1^4-2y1^3 y2+2y1 y2^3-y2^2((x1-x2)^2+y2^2)+x1 Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]-x2 Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]),{k1,k2},i},{((x1+x2)((x1-x2)^2+(y1-y2)^2)+Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)])/(2((x1-x2)^2+(y1-y2)^2)),1/(2((x1-x2)^2+(y1-y2)^2)(y1-y2))((x1-x2)^2 y1^2+y1^4-2y1^3 y2+2y1 y2^3-y2^2((x1-x2)^2+y2^2)-x1 Abs[y1-y2]Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]+Abs[y1-y2]x2 Sqrt[-((x1-x2)^2+(y1-y2)^2)((x1-x2)^2+(y1-y2)^2-4r^2)]),{k1,k2},i}},Nothing];
复制代码
定义两点连线的斜率函数s
  1. s[{x1_, y1_}, {x2_, y2_}] = (y1 - y2)/(x1 - x2);
复制代码
将所有交点存到intersec
  1. intersec = Select[Flatten[Table[(f /@ pair) /. r -> 2 Cos[\[Pi]/(2n) + i \[Pi]/n], {i, 0,(n-3)/2}], 2],! MemberQ[Table[MapThread[Chop@*N@*Plus, {-#[[1;;2]], Delete[points[[m]], 3]}], {m,n}], {0, 0}] &];
复制代码
把x轴上的交点存到ve
  1. ve = Select[
  2.   Flatten[Table[Thread[{(x /. Solve[(x - Cos[2 Pi i/n])^2 + (Sin[2 Pi i/n])^2 == (2 Cos[\[Pi]/(2n) + j \[Pi]/n])^2 && x > 0, x, Reals]), i, j}], {i, 0, (n-1)/2}, {j, 0, (n-3)/2}], 2],
  3.   Length[#] > 1 &];
复制代码
把x轴上方的交点存到nve
  1. nve = {}; For[i = 1, i < Length[intersec], i++,
  2. If[N[intersec[[i]][[2]]] > 0, AppendTo[nve, intersec[[i]]], ,
  3. Print[i]]]
复制代码
把倾斜角存到angle
  1. angle = Flatten[
  2.   Table[If[Chop[N[First[ve[[i]]] - First[nve[[j]]]]] !=
  3.      0, {Chop[N[ArcTan[s[{First[ve[[i]]], 0}, nve[[j]][[1 ;; 2]]]]]],
  4.      Flatten[{Delete[ve[[i]], 1], nve[[j]][[3 ;; 4]]}]}, Nothing], {i,
  5.      Length[ve]}, {j, Length[nve]}], 1];
复制代码
开始搜索
  1. S = {}; For[i = 1, i != Length[angle], i++,
  2. For[j = i + 1, j != Length[angle], j++,
  3.   If[Chop[N[
  4.       10^(-2) FractionalPart[(angle[[i]] + angle[[j]]) 2 n/Pi]]] == 0,
  5.     AppendTo[
  6.     S, {Last[angle[[i]]], Last[angle[[j]]],
  7.      Chop[N[(angle[[i]] + angle[[j]]) 2 n/Pi]]}]]]]
复制代码
加一个进度条
  1. Dynamic[{i, j}]
复制代码
加一个实时输出栏
  1. Dynamic[S]
复制代码
把+改成-再搜索一波
  1. S = {}; For[i = 1, i != Length[angle], i++,For[j = i + 1, j != Length[angle], j++, If[Chop[N[10^(-2) FractionalPart[(angle[[i]] - angle[[j]]) 2 n/Pi]]] == 0, AppendTo[S, {Last[angle[[i]]], Last[angle[[j]]], Chop[N[(angle[[i]] - angle[[j]]) 2 n/Pi]]}]]]]
复制代码
分享到: QQ空间QQ空间 腾讯微博腾讯微博 腾讯朋友腾讯朋友

测试一下代码
  1. ve[[50 ;; 60]]//N
复制代码
输出{{2.65636, 3., 0.}, {2.62782, 3., 1.}, {2.57098, 3., 2.}, {2.48631,
  3., 3.}, {2.37441, 3., 4.}, {2.23588, 3., 5.}, {2.07108, 3.,
  6.}, {1.87939, 3., 7.}, {1.65748, 3., 8.}, {0.139541, 3.,
  9.}, {1.39255, 3., 9.}}
  1. nve[[50 ;; 60]] // N
复制代码
输出{{2.49516, 1.6411, {2., 3., 0.}}, {2.26466,
  1.90027, {2., 4., 0.}}, {1.99477, 2.11433, {2., 5., 0.}}, {1.69543,
  2.27735, {2., 6., 0.}}, {1.3775, 2.3859, {2., 7., 0.}}, {1.05213,
  2.43911, {2., 8., 0.}}, {0.730101,
  2.43871, {2., 9., 0.}}, {0.421195,
  2.38872, {2., 10., 0.}}, {0.133676,
  2.29513, {2., 11., 0.}}, {-0.126118,
  2.16536, {2., 12., 0.}}, {-0.353996, 2.00761, {2., 13., 0.}}}
开始搜索(预计时长105天,我用3小时测试了一下,只搜索了大概千分之六,没搜索到任何结果)

TOP

回复 2# hbghlyj
虽然程序是正确的,但是n=15和n=9都没有搜索到任何结果,还是弃坑了.
这样搜索也是太投机主义了(捂脸

TOP

返回列表 回复 发帖