按弦长相等平分椭圆应该怎么做？

Question

比如：按弦长相等7等分椭圆

一种办法是以每段弦长所对应的圆心角为未知数建立联立方程组求解！

代码：

ClearAll["Global`*"]
a0 = 4;
b0 = 3;
r[t_] = Solve[1/a0^2 r^2 Cos[t]^2 + 1/b0^2 r^2 Sin[t]^2 == 1, r][[2, 
   1, 2]];

chord[i_] := Sqrt[
 r[Sum[a[k], {k, 1, i}]]^2 + r[Sum[a[k], {k, 1, i - 1}]]^2 - 
  2 r[Sum[a[k], {k, 1, i}]] r[Sum[a[k], {k, 1, i - 1}]] Cos[a[i]]]
divide[n_] := 
 FindRoot[{Equal @@ # & /@ Partition[Table[chord[i], {i, n}], 2, 1], 
   Sum[a[k], {k, 1, n}] == 2.0 Pi}, 
  Prepend[Table[{a[i], 2. Pi/n}, {i, 2, n}], {a[1], 1.0 Pi/n}]]
plot[n_] := 
 Show[ListLinePlot[(r[#] AngleVector[#] & /@ 
      Accumulate@divide[n][[All, 2]]) /. {x_, y__} -> {x, y, x}], 
  ParametricPlot[{a0 Cos[t], b0 Sin[t]}, {t, 0, 2 Pi}], 
  PlotRange -> All, AspectRatio -> Automatic]

运行结果：

{plot[3], plot[4], plot[5]}

但当要平分的段数太多时，这种办法效率非常低！

比如：

plot[400] // AbsoluteTiming

运行了十几分钟，除了这种方法，还有什么高效率的方法吗？

 

Answer 1

抱歉之前没看清题目，弦长相等的话就只能FindRoot了。

优化主要体现在两方面，一是利用椭圆的轴对称性分奇偶处理只计算上半圆，二是利用FindRoot的向量模式整体求解。除此之外就是一些比较细节的地方，比如去掉Sqrt之类。

代码如下，400用时1秒多一点。进一步优化的话可以加上编译，不过感觉这个效率已经可以接受了（编译过的代码贴在下面了，400用时0.05s）

eqnEvenQ[\[Theta]list_List] := Module[{lenList},
   lenList = 
    2 (a0^2 + b0^2 + (-a0^2 + b0^2) Cos[#1 + #2]) Sin[(#1 - #2)/
        2]^2 & @@@ Partition[N@Flatten@{0, \[Theta]list, \[Pi]}, 2, 1];
   lenList[[2 ;; -1]] - lenList[[1 ;; -2]]
   ];
eqnOddQ[\[Theta]list_List] := Module[{lenList},
   lenList = 
    2 (a0^2 + b0^2 + (-a0^2 + b0^2) Cos[#1 + #2]) Sin[(#1 - #2)/
        2]^2 & @@@ Partition[N@Flatten@{0, \[Theta]list}, 2, 1];
   AppendTo[lenList, (2 b0 Sin[\[Theta]list[[-1]]])^2];
   lenList[[2 ;; -1]] - lenList[[1 ;; -2]]
   ];
plot[n_?EvenQ] := Module[{\[Theta]list},
   \[Theta]list = 
    Flatten[{0, \[Theta] /. 
       FindRoot[
        eqnEvenQ[\[Theta]], {\[Theta], 
         Range[n/2 - 1]/(n/2)*\[Pi]}], \[Pi]}];
   ParametricPlot[{a0 Cos[t], b0 Sin[t]}, {t, 0, 2 Pi}, 
    Epilog -> {Red, 
      Line[{a0 Cos[#], b0 Sin[#]} & /@ 
        Join[\[Theta]list, -Reverse@\[Theta]list[[1 ;; -2]]]]}, 
    AspectRatio -> Automatic]
   ];
plot[n_?OddQ] := Module[{\[Theta]list},
   \[Theta]list = 
    Flatten[{0, \[Theta] /. 
       FindRoot[
        eqnOddQ[\[Theta]], {\[Theta], 
         Range[(n - 1)/2]/((n - 1)/2)*\[Pi]}]}];
   ParametricPlot[{a0 Cos[t], b0 Sin[t]}, {t, 0, 2 Pi}, 
    Epilog -> {Red, 
      Line[{a0 Cos[#], b0 Sin[#]} & /@ 
        Join[\[Theta]list, -Reverse@\[Theta]list]]}, 
    AspectRatio -> Automatic]
   ];
plot[400] // AbsoluteTiming

把编译过的也贴上来吧，400用时0.05s

eqnEvenQ[\[Theta]list_List] := cEven[\[Theta]list];
cEven = With[{a0 = a0, b0 = b0},
   Compile[{{\[Theta], _Real, 1}},
    Module[{lenList, \[Theta]list},
     \[Theta]list = Append[Prepend[\[Theta], 0], \[Pi]];
     lenList = 
      Table[2 (a0^2 + 
          b0^2 + (-a0^2 + 
             b0^2) Cos[\[Theta]list[[i]] + \[Theta]list[[
              i + 1]]]) Sin[(\[Theta]list[[i]] - \[Theta]list[[
           i + 1]])/2]^2, {i, Length[\[Theta]list] - 1}];
     lenList[[2 ;; -1]] - lenList[[1 ;; -2]]
     ], CompilationTarget -> "C"]
   ];
eqnOddQ[\[Theta]list_List] := cOdd[\[Theta]list];
cOdd = With[{a0 = a0, b0 = b0},
   Compile[{{\[Theta], _Real, 1}},
    Module[{lenList, \[Theta]list},
     \[Theta]list = Prepend[\[Theta], 0];
     lenList = 
      Table[2 (a0^2 + 
          b0^2 + (-a0^2 + 
             b0^2) Cos[\[Theta]list[[i]] + \[Theta]list[[
              i + 1]]]) Sin[(\[Theta]list[[i]] - \[Theta]list[[
           i + 1]])/2]^2, {i, Length[\[Theta]list] - 1}];
     AppendTo[lenList, (2 b0 Sin[\[Theta]list[[-1]]])^2];
     lenList[[2 ;; -1]] - lenList[[1 ;; -2]]
     ], CompilationTarget -> "C"]
   ];
plot[n_?EvenQ] := Module[{\[Theta]list},
   \[Theta]list = 
    Flatten[{0, \[Theta] /. 
       FindRoot[
        eqnEvenQ[\[Theta]], {\[Theta], 
         Range[n/2 - 1]/(n/2)*\[Pi]}], \[Pi]}];
   ParametricPlot[{a0 Cos[t], b0 Sin[t]}, {t, 0, 2 Pi}, 
    Epilog -> {Red, 
      Line[{a0 Cos[#], b0 Sin[#]} & /@ 
        Join[\[Theta]list, -Reverse@\[Theta]list[[1 ;; -2]]]]}, 
    AspectRatio -> Automatic]
   ];
plot[n_?OddQ] := Module[{\[Theta]list},
   \[Theta]list = 
    Flatten[{0, \[Theta] /. 
       FindRoot[
        eqnOddQ[\[Theta]], {\[Theta], 
         Range[(n - 1)/2]/((n - 1)/2)*\[Pi]}]}];
   ParametricPlot[{a0 Cos[t], b0 Sin[t]}, {t, 0, 2 Pi}, 
    Epilog -> {Red, 
      Line[{a0 Cos[#], b0 Sin[#]} & /@ 
        Join[\[Theta]list, -Reverse@\[Theta]list]]}, 
    AspectRatio -> Automatic]
   ];
plot[400] // AbsoluteTiming

==============================以下原答案=============================

可以用椭圆积分EllipticE，然后利用InverseFunction代替FindRoot。函数的具体说明参考帮助文档，我这里就只发代码了

a0 = 4;
b0 = 3;
plot[n_] := Module[{\[Theta]list},
\[Theta]list = 
InverseFunction[
EllipticE[#, 1 - a0^2/b0^2] &] /@ (N@
EllipticE[2 \[Pi], 1 - a0^2/b0^2]*Range[n]/n);
ParametricPlot[{a0 Cos[t], b0 Sin[t]}, {t, 0, 2 Pi}, 
Epilog -> {Red, Line[{a0 Cos[#], b0 Sin[#]} & /@ \[Theta]list]}, 
AspectRatio -> Automatic]
];
plot[400] // AbsoluteTiming

400份在我这里用时0.3秒左右

Comments

可不可以把方程参数化，以弦长为t

---

为什么一定要以弦长为参数？这种做法有哪里不满足要求吗？

---

可以的话，应该可以简单许多

---

想了下没太想到该怎么用弦长。。。这个方法显得麻烦主要是因为分了奇偶两部分，如果直接合到一起的话就比较短了

---

FindRoot中的变量t[1], t[2], t[3]...可以编译吗？（ps: 回复中的代码能高亮吗？）
a0 = 4.0;
b0 = 3.0;
t[0] = 0.0;
r2[t1_] = b0^2/(a0 + Sqrt[a0^2 - b0^2] Cos[t1]);
eqns[n_] :=
 Append[Apply[Equal,
   Partition[
    Table[r2[t[i]]^2 + r2[t[i - 1]]^2 -
      2 r2[t[i]] r2[t[i - 1]] Cos[t[i] - t[i - 1]], {i, n}], 2,
    1], {1}], t[1] + t[n - 1] == 2 \[Pi]]
plot2[n_] :=
 ListPlot[({Sqrt[a0^2 - b0^2] + r2[#1] Cos[#1], r2[#1] Sin[#1]} &) /@
   FindRoot[eqns[n], Evaluate[Table[{t[i], (i 2 \[Pi])/n}, {i, n}]]][[
    All, 2]], Joined -> True, AspectRatio -> Automatic]
plot2[400] // AbsoluteTiming