公告:1)网站程序升级:Q2A升级到1.8.6,Wordpress升级到5.7.2
2)修复了头像加载慢与提交问题反应慢等问题
2021-06-16

欢迎来到 Mathematica 问答社区

提问时请贴上文本代码

语法高亮:在编辑器中点击

被禁止的话题:广告破解

请阅读:《提问的智慧》

备用域名:mma.ooo

支持LaTex数学公式
行内公式标识符:\$ 或“$\backslash ($”+“$\backslash )$”,
行间公式标识符:\$\$ 或 “$\backslash [$”+“$\backslash ]$”

社区建议QQ群:365716997

分类

+2 投票
1.2k 浏览

比如:按弦长相等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

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

 

分类:绘图 | 用户: niturpe (251 分)
重新分类 用户:niturpe

1个回答

+4 投票
 
已采纳

抱歉之前没看清题目,弦长相等的话就只能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秒左右

用户: 无影东瓜 (1.2k 分)
采纳于 用户:niturpe
可不可以把方程参数化,以弦长为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
...