数学中国

标题: 椭圆周长数值积分与近似公式计算精度比较 [打印本页]

作者: Ysu2008 时间: 2023-3-20 19:53
标题: 椭圆周长数值积分与近似公式计算精度比较
本帖最后由 Ysu2008 于 2023-3-20 19:59 编辑

项名达级数公式\(C_0=2\pi a\left[ 1-\left( \frac{1}{2}\right)^2e^2-\left( \frac{1\cdot3}{2\cdot4}\right)^2\frac{e^4}{3}-\left( \frac{1\cdot3\cdot5}{2\cdot4\cdot6}\right)^2\frac{e^6}{5}-\left( \frac{1\cdot3\cdot5\cdot7}{2\cdot4\cdot6\cdot8}\right)^2\frac{e^8}{7}-......\right]\)
取项名达级数公式前100万项截断，作为“精确值”，计为\(C_0\)；

龙贝格数值积分此式 \(C_1=4\int_0^{\frac{\pi}{2}}\sqrt{a^2\sin^2t+b^2\cos^2t}dt\)

周长近似公式为 \(\begin{align}
C_2 &\approx \pi \left( {a + b} \right)\left[ {1 + \frac{{3{\lambda ^2}}}{{10 + \sqrt {4 - 3{\lambda ^2}} }} + \frac{3}{{{2^{17}}}}(1 + \frac{{5767168 - 1835041\pi }}{{33\pi }}{\lambda ^{4.88 + 11.04{\lambda ^{4.74}}}}){\lambda ^{10}}} \right]
\end{align}\)

取椭圆长半轴固定为\(a=1\)，短半轴从\(0.01\to1\)，步进\(0.01\)，总计100个椭圆周长，比较\(\left| C_0-C_1\right|\)与\(\left| C_0-C_2\right|\)精度大小。

[数据]：

编号	长半轴	短半轴	离心率e	项名达周长	数值积分周长	近似公式周长	数值积分误差	近似公式误差
1	1	0.01	0.9999	4.00109833	4.00109833	4.001144916	1.17E-12	4.66E-05
2	1	0.02	0.9998	4.00383916	4.00383916	4.003866361	6.13E-13	2.72E-05
3	1	0.03	0.9995	4.00790945	4.00790945	4.007921397	1.22E-13	1.19E-05
4	1	0.04	0.9992	4.013143313	4.013143313	4.013147132	2.02E-13	3.82E-06
5	1	0.05	0.9987	4.019425619	4.019425619	4.019425851	1.30E-13	2.31E-07
6	1	0.06	0.9982	4.026668233	4.026668233	4.026667256	4.99E-13	9.78E-07
7	1	0.07	0.9975	4.034799933	4.034799933	4.03479884	5.05E-13	1.09E-06
8	1	0.08	0.9968	4.043761127	4.043761127	4.043760335	1.10E-13	7.92E-07
9	1	0.09	0.9959	4.053500734	4.053500734	4.053500336	2.55E-13	3.98E-07
10	1	0.1	0.9950	4.06397418	4.06397418	4.063974129	4.62E-14	5.13E-08
11	1	0.11	0.9939	4.075142024	4.075142024	4.075142229	1.01E-13	2.04E-07
12	1	0.12	0.9928	4.086968989	4.086968989	4.086969356	4.44E-15	3.67E-07
13	1	0.13	0.9915	4.09942324	4.09942324	4.099423691	2.22E-14	4.52E-07
14	1	0.14	0.9902	4.112475833	4.112475833	4.112476311	1.49E-13	4.78E-07
15	1	0.15	0.9887	4.126100294	4.126100294	4.126100759	2.14E-13	4.65E-07
16	1	0.16	0.9871	4.140272276	4.140272276	4.140272703	2.33E-13	4.27E-07
17	1	0.17	0.9854	4.154969285	4.154969285	4.154969662	1.03E-13	3.77E-07
18	1	0.18	0.9837	4.170170459	4.170170459	4.170170781	2.54E-13	3.22E-07
19	1	0.19	0.9818	4.185856378	4.185856378	4.185856645	1.21E-13	2.67E-07
20	1	0.2	0.9798	4.202008908	4.202008908	4.202009124	2.61E-13	2.16E-07
21	1	0.21	0.9777	4.21861107	4.21861107	4.218611241	4.99E-13	1.71E-07
22	1	0.22	0.9755	4.235646924	4.235646924	4.235647057	2.31E-14	1.33E-07
23	1	0.23	0.9732	4.253101475	4.253101475	4.253101575	7.82E-14	1.01E-07
24	1	0.24	0.9708	4.270960581	4.270960581	4.270960656	1.31E-13	7.45E-08
25	1	0.25	0.9682	4.289210888	4.289210888	4.289210941	5.42E-14	5.37E-08
26	1	0.26	0.9656	4.307839754	4.307839754	4.307839791	5.33E-14	3.76E-08
27	1	0.27	0.9629	4.3268352	4.3268352	4.326835225	3.53E-13	2.52E-08
28	1	0.28	0.9600	4.346185854	4.346185854	4.34618587	6.84E-14	1.60E-08
29	1	0.29	0.9570	4.365880906	4.365880906	4.365880916	2.18E-13	9.36E-09
30	1	0.3	0.9539	4.38591007	4.38591007	4.385910074	3.92E-13	4.64E-09
31	1	0.31	0.9507	4.406263539	4.406263539	4.406263541	7.82E-14	1.42E-09
32	1	0.32	0.9474	4.426931963	4.426931963	4.426931962	2.72E-13	6.71E-10
33	1	0.33	0.9440	4.447906407	4.447906407	4.447906405	6.22E-14	1.94E-09
34	1	0.34	0.9404	4.469178334	4.469178334	4.469178332	3.38E-13	2.63E-09
35	1	0.35	0.9367	4.490739573	4.490739573	4.49073957	3.55E-14	2.91E-09
36	1	0.36	0.9330	4.512582298	4.512582298	4.512582295	4.00E-14	2.93E-09
37	1	0.37	0.9290	4.534699008	4.534699008	4.534699005	3.97E-13	2.77E-09
38	1	0.38	0.9250	4.557082509	4.557082509	4.557082506	8.88E-16	2.52E-09
39	1	0.39	0.9208	4.579725891	4.579725891	4.579725889	2.49E-14	2.23E-09
40	1	0.4	0.9165	4.602622519	4.602622519	4.602622517	2.06E-13	1.93E-09
41	1	0.41	0.9121	4.625766012	4.625766012	4.625766011	3.55E-13	1.63E-09
42	1	0.42	0.9075	4.649150232	4.649150232	4.649150231	3.55E-13	1.36E-09
43	1	0.43	0.9028	4.672769271	4.672769271	4.67276927	7.99E-15	1.12E-09
44	1	0.44	0.8980	4.696617437	4.696617437	4.696617436	2.02E-13	9.08E-10
45	1	0.45	0.8930	4.720689244	4.720689244	4.720689243	9.24E-14	7.29E-10
46	1	0.46	0.8879	4.744979404	4.744979404	4.744979404	7.99E-15	5.79E-10
47	1	0.47	0.8827	4.769482814	4.769482814	4.769482813	7.02E-14	4.56E-10
48	1	0.48	0.8773	4.794194546	4.794194546	4.794194546	1.59E-13	3.55E-10
49	1	0.49	0.8717	4.819109844	4.819109844	4.819109844	2.40E-13	2.75E-10
50	1	0.5	0.8660	4.84422411	4.84422411	4.84422411	2.66E-14	2.11E-10
51	1	0.51	0.8602	4.869532901	4.869532901	4.8695329	4.49E-13	1.60E-10
52	1	0.52	0.8542	4.895031917	4.895031917	4.895031917	3.85E-13	1.21E-10
53	1	0.53	0.8480	4.920717002	4.920717002	4.920717001	9.15E-14	9.04E-11
54	1	0.54	0.8417	4.946584128	4.946584128	4.946584128	2.53E-13	6.71E-11
55	1	0.55	0.8352	4.9726294	4.9726294	4.9726294	3.41E-13	4.94E-11
56	1	0.56	0.8285	4.998849042	4.998849042	4.998849042	1.69E-14	3.61E-11
57	1	0.57	0.8216	5.025239393	5.025239393	5.025239393	2.50E-13	2.62E-11
58	1	0.58	0.8146	5.051796909	5.051796909	5.051796909	2.40E-13	1.88E-11
59	1	0.59	0.8074	5.078518147	5.078518147	5.078518147	1.24E-13	1.34E-11
60	1	0.6	0.8000	5.105399773	5.105399773	5.105399773	2.66E-13	9.51E-12
61	1	0.61	0.7924	5.132438546	5.132438546	5.132438546	3.72E-13	6.67E-12
62	1	0.62	0.7846	5.159631325	5.159631325	5.159631325	3.02E-14	4.63E-12
63	1	0.63	0.7766	5.186975056	5.186975056	5.186975056	2.86E-13	3.19E-12
64	1	0.64	0.7684	5.214466775	5.214466775	5.214466775	3.23E-13	2.18E-12
65	1	0.65	0.7599	5.2421036	5.2421036	5.2421036	2.65E-13	1.47E-12
66	1	0.66	0.7513	5.269882731	5.269882731	5.269882731	2.19E-13	9.80E-13
67	1	0.67	0.7424	5.297801446	5.297801446	5.297801446	4.17E-14	6.48E-13
68	1	0.68	0.7332	5.325857098	5.325857098	5.325857098	2.66E-13	4.24E-13
69	1	0.69	0.7238	5.354047112	5.354047112	5.354047112	1.81E-13	2.74E-13
70	1	0.7	0.7141	5.382368981	5.382368981	5.382368981	3.11E-14	1.74E-13
71	1	0.71	0.7042	5.410820269	5.410820269	5.410820269	8.88E-15	1.10E-13
72	1	0.72	0.6940	5.4393986	5.4393986	5.4393986	5.32E-13	6.93E-14
73	1	0.73	0.6834	5.468101664	5.468101664	5.468101664	2.66E-13	4.09E-14
74	1	0.74	0.6726	5.496927208	5.496927208	5.496927208	2.49E-14	2.58E-14
75	1	0.75	0.6614	5.52587304	5.52587304	5.52587304	3.36E-13	1.51E-14
76	1	0.76	0.6499	5.554937023	5.554937023	5.554937023	2.03E-13	7.99E-15
77	1	0.77	0.6380	5.584117073	5.584117073	5.584117073	1.06E-13	5.33E-15
78	1	0.78	0.6258	5.613411161	5.613411161	5.613411161	2.74E-13	1.78E-15
79	1	0.79	0.6131	5.642817306	5.642817306	5.642817306	9.59E-14	1.78E-15
80	1	0.8	0.6000	5.672333578	5.672333578	5.672333578	3.73E-14	8.88E-16
81	1	0.81	0.5864	5.701958092	5.701958092	5.701958092	4.99E-13	8.88E-16
82	1	0.82	0.5724	5.731689012	5.731689012	5.731689012	1.27E-13	8.88E-16
83	1	0.83	0.5578	5.761524543	5.761524543	5.761524543	1.04E-13	0.00E+00
84	1	0.84	0.5426	5.791462935	5.791462935	5.791462935	3.02E-14	8.88E-16
85	1	0.85	0.5268	5.82150248	5.82150248	5.82150248	2.93E-13	0.00E+00
86	1	0.86	0.5103	5.851641509	5.851641509	5.851641509	8.97E-14	8.88E-16
87	1	0.87	0.4931	5.881878392	5.881878392	5.881878392	4.10E-13	8.88E-16
88	1	0.88	0.4750	5.912211538	5.912211538	5.912211538	2.63E-13	8.88E-16
89	1	0.89	0.4560	5.942639391	5.942639391	5.942639391	1.35E-13	0.00E+00
90	1	0.9	0.4359	5.973160433	5.973160433	5.973160433	2.34E-13	0.00E+00
91	1	0.91	0.4146	6.003773177	6.003773177	6.003773177	2.70E-13	8.88E-16
92	1	0.92	0.3919	6.034476173	6.034476173	6.034476173	4.72E-13	8.88E-16
93	1	0.93	0.3676	6.065268001	6.065268001	6.065268001	1.60E-13	8.88E-16
94	1	0.94	0.3412	6.096147275	6.096147275	6.096147275	5.34E-13	0.00E+00
95	1	0.95	0.3122	6.127112637	6.127112637	6.127112637	3.71E-13	1.78E-15
96	1	0.96	0.2800	6.158162759	6.158162759	6.158162759	2.18E-13	0.00E+00
97	1	0.97	0.2431	6.189296344	6.189296344	6.189296344	1.71E-13	0.00E+00
98	1	0.98	0.1990	6.220512121	6.220512121	6.220512121	4.41E-13	8.88E-16
99	1	0.99	0.1411	6.251808848	6.251808848	6.251808848	1.03E-12	8.88E-16
100	1	1	0.0000	6.283185307	6.283185307	6.283185307	8.88E-16	0.00E+00

[结论]
[attach]124977[/attach]
椭圆离心率越大（越接近1），近似公式误差越大。大约在\(e\approx0.75\)时，近似公式与龙贝格积分持平，越接近0精度越高。

数值积分的优势在于各种离心率下表现都比较稳定，而且可以计算指定区域的椭圆弧长。

作者: uk702 时间: 2023-3-20 20:16
本帖最后由 uk702 于 2023-3-28 12:56 编辑

xxxxxxxxxxxxxxxxxxxxx

作者: uk702 时间: 2023-3-20 21:22
本帖最后由 uk702 于 2023-3-28 12:56 编辑

xxxxxxxxxxxxxxxxxx

作者: elim 时间: 2023-3-20 23:10
主贴的结果非常漂亮, 谢谢 Ysu2008 !

从这些数值计算结果可以轻易弄出一个优于简单而已有结果的近似公式来。关于这点，各位有什么怀疑吗？

作者: elim 时间: 2023-3-21 00:37

永远发表于 2023-3-20 08:21
各位帮我看看e老师这个拉马努金公式补尝拟合函数族是怎么构造出来的，谁会？？？

这个函数族(u,v 视为参量）是我提出来的。但我没有说过这是拉氏拟合的补偿拟合。

作者: Ysu2008 时间: 2023-3-21 09:33

uk702 发表于 2023-3-20 21:22
若楼主有时间，望楼主帮忙检查核对一下，下面的“简单”公式，最大误差不超过 5.54*10^-6 。

\(a=1{,}b=0.0034\)
项名达公式100万截断周长\(C_0\ =4.000151905348548\)
你的公式算得周长 \(C_2=4.0001574446194175\)
误差\(\left| C_0-C_2\right|=5.539270869547863E-6\)

作者: uk702 时间: 2023-3-21 09:57
本帖最后由 uk702 于 2023-3-28 12:56 编辑

xxxxxxxxxxxxxxxxxxx

作者: dodonaomikiki 时间: 2023-3-21 15:00
zhichi支持！
非常漂亮！
赞美

作者: elim 时间: 2023-3-23 09:48
我们所考虑过的补偿拟合的误差，差不多在 \(10^{-5}\) 的数量级. 这是有原因的。
椭圆周长的多项式拟合的极致就是超几何级数。使用了无穷多个自由度(级数系数)
这就解释了寥寥几个参变量是不足以突破 \(10^{-5}\) 这个坎的。要突破，考虑以
下简单方案：
给定\(k\),取最小\(m\)使\(\;r_m\small=\displaystyle\sum_{n = m+1}^\infty\binom{1/2}{n}^2\le 10^{-k},\)令\(\small G(x)=\displaystyle\sum_{n=0}^m\binom{1/2}{n}^2x^{2n}\)
考虑 \(\small\dfrac{{\scriptsize\displaystyle\sum_{n = m+1}^\infty\binom{1/2}{n}^2}x^{2n}}{\binom{1/2}{m+1}^2 x^{2m+2}}\) 的形如 \(\small 1+\big(\dfrac{4}{\pi}-r_m\big)\varphi(x)\) 的拟合,其中\(\small0\le\varphi\le 1.\)

注意以下椭圆周长的高精度极速算法。数值积分望尘莫及
F(x)=my(t=x^2/4,s=1+t,k=1);if(x==1,return(4/Pi));while((t=t*x^2*((2*k-1)/(2*(k+1)))^2)&&(t> 10^(-100))&&(s=s+t),k=k+1);return(s);

作者: elim 时间: 2023-3-23 17:12
\(\small C=\pi(a+b)\big(1+\dfrac{3\lambda^2}{10+\sqrt{4-3\lambda^2}}+\dfrac{3}{2^{17}}\lambda^{10}\big(1+\dfrac{\mu\lambda}{(1+(1-\lambda^u)^v)^w}\big)\big)\)
\(\lambda={\large\frac{a-b}{a+b}},\;\epsilon=\small\dfrac{2\sqrt{\lambda}}{1+\lambda},\;\lambda=\dfrac{2-\epsilon^2+2\sqrt{1-\epsilon^2}}{\epsilon^2}\)
\(\mu=21.381505539,\;u=2.02777,\;v=0.823239,\;w=4.84343\)

作者: dodonaomikiki 时间: 2023-3-23 17:32

太佩服太佩服啦！
佩服E老师，这些东西，我在学校里从木有听说过~~~~~~可谓，孤陋寡闻也！数学专业之精深知识也！

作者: Ysu2008 时间: 2023-3-23 17:35

elim 发表于 2023-3-23 17:12
\(\small C=\pi(a+b)\big(1+\dfrac{3\lambda^2}{10+\sqrt{4-3\lambda^2}}+\dfrac{3}{2^{17}}\lambda^{10}\b ...

弧长近似值可不可以推？

作者: elim 时间: 2023-3-23 17:56
本帖最后由 elim 于 2023-3-23 03:12 编辑

elim 发表于 2023-3-23 02:12
\(\small C=\pi(a+b)\big(1+\dfrac{3\lambda^2}{10+\sqrt{4-3\lambda^2}}+\dfrac{3}{2^{17}}\lambda^{10}\b ...

楼上的拟合误差一致地在 \(\le 4.3\times 10^{-6}\) 范围。三个参数 \(u,v,w\), 已无多少优化余地。我用了以下的梯度优化办法。

[attach]125103[/attach]

作者: elim 时间: 2023-3-23 21:01

永远发表于 2023-3-23 04:34
请问e老师，这两个参数有什么用处？？？

\(\epsilon=\small\dfrac{2\sqrt{\lambda}}{1+\lambda},\;\ ...

离心率与 \(\lambda\) 的互换。

作者: elim 时间: 2023-3-23 23:13
本帖最后由 elim 于 2023-3-23 14:56 编辑

我就是好奇推了一下 \(\lambda=(a-b)/(a+b),\;\epsilon = \sqrt{1-(b/a)^2}\) 之间的关系。你用不着不用就跟我急.

a, b, c 初始值需要先输入。
所有定义的函数需要依次执行一下。

作者: elim 时间: 2023-3-24 06:02

永远发表于 2023-3-23 08:21
又试了一下，不明所以然，还是不行，初始值怎么先输入，定义函数怎么依次执行？？？我22楼已上传好文件 ...

我自己的文件，每次打开都要激活定义：把鼠标置于定义，shift+return
初始值可以用目测得到的那些值。

建议使用 gp/pari 作这种运算，或请教天山草老师。

作者: elim 时间: 2023-3-24 15:26
本帖最后由 elim 于 2023-3-24 00:38 编辑

永远发表于 2023-3-23 15:55
请老师把你的带有初始值的mathematica文件发出来看看

[attach]125179[/attach]
这是我现在运行的nb 文件. 当初在 a = a -t u[a,b,c] 那一行就是 a = 1.8... 把这些赋值完成后，就运行 a = a -t u[a,b,c] ...
你现在的问题与椭圆无关，是怎么用 Mathematica 的问题。我不常用它。帮不了你.

作者: elim 时间: 2023-3-24 21:19
你的程序的错误是怎么产生，又是怎么搞定的？你的‘事实上’的事实是什么？弄懂我的不依赖目测确定参量的梯度法了吗？可否依样画葫芦弄个四参量的版本看看？
我的建议是 \(\small\sigma(x)=1+{\large\frac{(\mu-1) x^2}{(1+a(1-x^b)^c)^d}},\;\mu=\frac{2^{17}}{3}\big(\frac{4}{\large\pi}-\frac{14}{11}\big)\)

作者: elim 时间: 2023-3-24 21:32
你的能力的限度是什么？不能学习了？

作者: elim 时间: 2023-3-25 00:06
回答永远你的问题之前总有点纠结：回答有什么用，还会不会引出更多的无聊问题．
我来问个问题吧：兰登原文名是什么？他的变换公式的证明在哪里？

作者: elim 时间: 2023-3-25 06:53
问你兰登原文名是什么？他的变换公式的证明在哪里？

作者: elim 时间: 2023-3-25 21:07

永远发表于 2023-3-24 08:01
均方误差相关疑问

请问程序第一行p[n_]：=1-(1-n/20)^2其中“1-(1-n/20)^2”是怎么来的，这个数学表 ...

这不是均方误差问题。而是多点误差平方和的最小化问题。求出这种和相对于（u,v,w) 的梯度，优化（u,v,w).
使得这个差方和趋于最小。为了在拟合误差较大之处加权，就不采用均匀分布的计算点了。

作者: elim 时间: 2023-3-25 21:13
你不知道我在干什么，问我第一步要怎么干有意思吗？

作者: elim 时间: 2023-4-1 09:13
本贴初略探讨椭圆周长的拟合函数的一种构造方案和误差。
[attach]125563[/attach]

\(x(1-x)\) 是一个两端为\(0\),其他处为正的函数，\(x^u(1-x)^w\) 对\(u,w>0\)有相同的性质，但改变了两端的曲率甚至符号, 加上参数\(v\)进一步对\(x^u(1-x^v)^w\)的最大值位置有了控制，最后，参数\(z\)使得\((1+(x^u(1-x^v)^w)^z\)的峰值有明显的控制。

欢迎光临数学中国 (http://www.mathchina.com/bbs/)