用Mathematica编程求出最大误差函数值时如：10⁻⁸下的(x,a,b,c,d)

笨笨

Jack315 发表于 2023-9-14 09:51
拟合函数：\(\lambda_0+\lambda_1x+\lambda_2x^2\)
窗函数为：\(\frac{1}{2tan^{-1}(\frac{b-a}{2}k)}\bi ...

代码不是楼主专业，若代码有疑问，请把代码发出来，看看路过的高手能不能解决一下？

笨笨

elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨，而是刻意不动脑子不学习不动手的笨．

@Jack315 :  建议用切比雪夫理论一举得到支集 ...

elim老师自己说的方案，可能自己都实现不了？或者只是一种猜想。每次回帖都是留下一句话，然后就没有然后了，要是能给个具体实例分析一下就好了。可惜了，出于尊重，楼主持保留意见。

笨笨

事实上，增加有限项根本提高不了精度，下面就是一个失败的案例


\(\large\mu  = \frac{{{\rm{5905580032 - 1879201796}}\pi }}{{{\rm{64559}}\pi }}\)
\( F\left( x \right) - 1 - \frac{{3{x^2}}}{{10 + \sqrt {4 - 3{x^2}} }} - \frac{3}{{{2^{17}}}}{x^{10}} -\frac{{79}}{{{2^{21}}}}{x^{12}} -\frac{{1459}}{{{2^{25}}}}{x^{14}} - \frac{{5869}}{{{2^{27}}}}{x^{16}}\left[ {1 + \frac{{\left( {\frac{{{\rm{353055}}}}{{{\rm{375616}}}} + \left( {\mu  - 1 - \frac{{{\rm{353055}}}}{{{\rm{375616}}}}} \right)x} \right){x^2}}}{{{{\left( {1 + {x^{{\rm{ - 3}}{\rm{.98356}}}}{{\left( {1 - {x^{15.98099}}} \right)}^{{\rm{0}}{\rm{.937338}}}}} \right)}^{{\rm{1}}{\rm{.016278}}}}}}} \right]\)

误差函数族界最小值：\(1.635594 \ldots  \times {10^{ - 7}},x \to 0.9297993309 \ldots \)

Jack315

笨笨 发表于 2023-9-14 21:42
代码不是楼主专业，若代码有疑问，请把代码发出来，看看路过的高手能不能解决一下？ ...

敢问 LZ 是什么专业 。言归正传……

计算误差函数数据的代码：

1. 计算误差函数数据[误差函数_] := Module[
   
2.   (*=======================================================*)
   
3.   (*输入参数：*)
   
4.   (*误差函数[x]=目标函数[x]-\!\(
   
5. \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(n\)]\(
   
6. \(\*SubscriptBox[\(拟合函数\), \(i\)]\)[x,
   
7. \*SubscriptBox[\(拟合参数\), \(i\)]]\)\)。*)
   
8.   (*输出数据：
   
9.   {{
   
10.   极值间低点{Subscript[x, 1],Subscript[x, 2],...,Subscript[x, m]},
    
11.   所有极值点{{Subscript[x, 1],Subscript[y, 1]},{Subscript[x, 2],Subscript[
    
12.   y, 2]},...,{Subscript[x, n],Subscript[y, n]}},
    
13.   极值范围{min,max},极值数量 n,
    
14.   最大极值点{x,y},
    
15.   待拟合点{Subscript[x, fit],Subscript[y, fit]},
    
16.   待拟合点两侧低点{lef,right},
    
17.   加窗(T/F)
    
18.   }}
    
19.   *)
    
20.   (*-------------------------------------------------------*)
    
21.   {位置,
    
22.    所有极值点X, 所有极值点Y, 极值范围, 极值数量,
    
23.    最大极值点X, 最大极值点Y,
    
24.    高度, 极值间低点, 低点分类,
    
25.    待拟合点X, 待拟合点Y,
    
26.    待拟合点两侧低点, 加窗
    
27.    }  (*end local variables*),
    
28.   (*=======================================================*)
    
29.   
    
30.   (*=======================================================*)
    
31.   (*代码体*)
    
32.   (*=======================================================*)
    
33.   
    
34.   (*-------------------------------------------------------*)
    
35.   (*计算 误差函数[x] 的极值点。*)
    
36.   (*-------------------------------------------------------*)
    
37.   所有极值点X = x /. NSolve[{D[误差函数[x], x] == 0, 0 <= x <= 1}, x];
    
38.   (*确定极值间低点的高度：不能包含任何极值点。*)
    
39.   (*可能存在的问题：极值点比极值间低点离 X 轴更近*)
    
40.   高度 = 0.95 Min[Abs[误差函数[所有极值点X]]];
    
41.   位置 = Length[所有极值点X];
    
42.   (*将边界点 x=0 和 x=1 当作极值点处理*)
    
43.   If[位置 == 0 || 所有极值点X[[1]] != 0, PrependTo[所有极值点X, 0]];
    
44.   If[位置 == 0 || 所有极值点X[[位置]] != 1, AppendTo[所有极值点X, 1]];
    
45.   所有极值点Y = N[误差函数[所有极值点X]];
    
46.   (*-------------------------------------------------------*)
    
47.   (*极值点数据统计。*)
    
48.   极值范围 = {Min[Abs[所有极值点Y]], Max[Abs[所有极值点Y]]};
    
49.   极值数量 = Length[所有极值点X];
    
50.   (*-------------------------------------------------------*)
    
51.   (*计算 |极值| 最大的点。*)
    
52.   最大极值点Y = 极值范围[[2]];
    
53.   最大极值点X = Part[所有极值点X, Part[Position[Abs[所有极值点Y], 最大极值点Y], 1, 1]];
    
54.   (*-------------------------------------------------------*)
    
55.   (*确定待拟合点数据。*)
    
56.   待拟合点X = 最大极值点X;
    
57.   待拟合点Y = N[误差函数[待拟合点X]];
    
58.   (*=======================================================*)
    
59.   (*-------------------------------------------------------*)
    
60.   (*计算 |误差函数[x]| 在相邻极值点区间内的低点。*)
    
61.   (*-------------------------------------------------------*)
    
62.   (*计算方程 误差函数[x] \[Equal] 0\[PlusMinus]高度 的解。*)
    
63.   极值间低点 = NSolve[{误差函数[x]^2 == 高度^2, 0 <= x <= 1}, x];
    
64.   If[Length[极值间低点] == 0,
    
65.    (*极值间不含任何低点。*)
    
66.    极值间低点 = {};
    
67.    待拟合点两侧低点 = {0, 1};
    
68.    加窗 = False,
    
69.    (*-------------------------------------------------------*)
    
70.    (*极值间至少含有一个低点。*)
    
71.    (*以极值点的 x 坐标做为分界线将 极值间低点 分类。*)
    
72.    (*分类结果：{"极值间低点"索引，"极值间低点"分类号}*)
    
73.    极值间低点 = x /. 极值间低点;
    
74.    低点分类 = Transpose[Table[
    
75.       Catch[
    
76.        For[j = 1, j < Length[所有极值点X], j++,
    
77.         If[所有极值点X[[j]] < 极值间低点[[i]] < 所有极值点X[[j + 1]], Throw[{i, j}]]
    
78.         ] (*end For*)
    
79.        ] (*end Catch*),
    
80.       {i, Length[极值间低点]}] (*end Table*)
    
81.      ];(*end Transpose*)
    
82.    (*-------------------------------------------------------*)
    
83.    (*根据低点分类，将极值间低点进行合并。*)
    
84.    极值间低点 = Table[
    
85.      Mean[
    
86.       Part[极值间低点, Flatten[Position[低点分类[[2]], i] (*end Position*)
    
87.         ] (*end Flatten*)
    
88.        ] (*end Part*)
    
89.       ] (*end Mean*),
    
90.      {i, Max[Part[低点分类, 2]]}] (*end Table*);
    
91.    (*-------------------------------------------------------*)
    
92.    (*极值间低点 在边界处修正。*)
    
93.    位置 = Length[极值间低点];
    
94.    If[Abs[误差函数[0]] < Abs[误差函数[极值间低点[[1]]]], 极值间低点[[1]] = 0];
    
95.    If[Abs[误差函数[1]] < Abs[误差函数[极值间低点[[位置]]]], 极值间低点[[位置]] = 1];
    
96.    (*-------------------------------------------------------*)
    
97.    (*确定窗函数数据。*)
    
98.    (*-------------------------------------------------------*)
    
99.    位置 = Part[Position[极值间低点, Part[
    
100.        Select[极值间低点, # > 待拟合点X &, 1], 1]], 1, 1];
     
101.    待拟合点两侧低点 = 极值间低点[[位置 - 1 ;; 位置]];
     
102.    加窗 = ! (待拟合点X == 0 || 待拟合点X == 1 || 待拟合点两侧低点 == {0, 1});
     
103.    ]; (*end If*)
     
104.   (*-------------------------------------------------------*)
     
105.   (*输出 误差函数[x] 数据。*)
     
106.   (*-------------------------------------------------------*)
     
107.   {{极值间低点,
     
108.     Transpose[{所有极值点X, 所有极值点Y}],
     
109.     极值范围, 极值数量,
     
110.     {最大极值点X, 最大极值点Y},
     
111.     {待拟合点X, 待拟合点Y},
     
112.     待拟合点两侧低点, 加窗
     
113.     }}
     
114.   ] (*end Module*)

复制代码

代码通过了下列两种情况的测试：
1. 误差函数为 2F1；
2. 误差函数为经部分拟合后的剩余误差。

误差函数经多次迭代后，极值点及其间低点究竟会是什么情况，
有点没把握。光靠想像很难判明情况。
如果上述两种情形不能完全覆盖所有情况，代码就可能会有问题。

各位路过的高手多提宝贵意见，谢过先

笨笨

Jack315 发表于 2023-9-15 15:09
敢问 LZ 是什么专业 。言归正传……

计算误差函数数据的代码：

感谢先生一路相伴，虽然得出的结果差强人意，但就凭这份执着的求真态度，令人敬佩，楼主目前人在外省，电脑不在身边，所以有啥想法也没法实现。

86楼是楼主的个人极限，希望志同道合者有人能突破这个极限。最好对拉马努金公式再拟合找到一个简洁的再拟合公式一举突破10^-12，那就牛逼了！

Jack315

笨笨 发表于 2023-8-30 17:20
这是楼主目前找的比较不错的参数与拟合函数族，小结一下：

\(\mu  = \left( {\frac{4}{\pi } - \frac{{14} ...

这是在 18# 前辈的基础上用 Sin 函数一次性再补偿的效果：





补偿函数为：
\[errMax\times\sin[\theta2x^{-1}(\theta)]\]

其中：
\(errMax\) 是 18# 前辈得到的结果：误差最大值： 1.82658e-7。
\(\theta2x(\theta)\) 是由误差函数的零点和极值点对应的角度转换为 x 值的函数：
\[\theta2x(\theta)=(0.643608 + 0.0137754 \theta - 0.0011511 \theta^2 + 0.0000270624 \theta^3) tan^{-1}[0.529581 (2.73772 + \theta)]\]

\(\theta2x(\theta)\)  函数有比较准确的关系，但实际需要的是其反函数。
用反函数表示的函数不知道能不能算。
由这个函数算出零度对应的位置为 ：\(x=0.622366\) 。
\(x=0\) 这个点是唯一不符合这个函数关系的异常点。
【猜想】这或许表示梯度迭代的天花板可能更低。

最后误差比较大的地方大约在 \(x<0.6\) 和 \(x>0.99\) 的区域。
最终仍然需要用加窗函数的办法来解决。

单峰函数和窗函数有比较多的细节问题……有时间再接着玩。

笨笨

Jack315 发表于 2023-9-20 16:01
这是在 18# 前辈的基础上用 Sin 函数一次性再补偿的效果：

看看86楼我的历史记录，能否突破？

笨笨

elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨，而是刻意不动脑子不学习不动手的笨．

@Jack315 :  建议用切比雪夫理论一举得到支集 ...

请问怎么用切比雪夫理论一举得到最佳多项式一致逼近？可否具体点？？？？

Jack315

笨笨 发表于 2023-9-20 20:14
看看86楼我的历史记录，能否突破？

稍微再改进了一点点：


最终误差公式：
\[\begin{align}
err(x)=&Hypergeometric2F1[-\frac{1}{2}, -\frac{1}{2}, 1, x^2]\notag\\&-(1+\frac{3 x^2}{10+\sqrt{4-3x^2}})\notag\\&-\frac{3 x^{10}}{131072} \bigg[1+\frac{\bigg(\frac{79}{48}+\big(-\frac{127}{48}+\frac{131072}{3}(-\frac{14}{11}+\frac{4}{\pi})\big) x\bigg)x^2}{\bigg(1+\frac{(1-x^{11.4687})^{0.937768}}{x^{2.6301}}\bigg)^{1.36059}}\bigg]\notag\\&-fit1(x)-fit2(x)\notag
\end{align}\]
其中：
\[\begin{align}
fit1(x)=&4.33118\times10^{-8} (0.917361+2.99355x-24.8171x^2+99.1052x^3\notag\\&-244.232x^4+396.284x^5-369.54x^6+137.966 x^7)\times\notag\\&[\tan^{-1}\big(1000000 (x+0.17365)\big)-\tan^{-1}\big(1000000(x-0.766799)\big)]\times\notag\\&[1-\sin(0.145098 - 3.65847 x - 4.21967 x^2 + 0.951311 x^3)]-\notag\\&2.72031\times10^{-7}\sin\big[InverseFunction\notag\\&\big[\tan^{-1}[0.529581 (2.73772 + x)]\times(0.643608 + 0.0137754x-0.0011511 x^2 + 0.0000270624 x^3)\big]\big]\notag
\end{align}\]

\[\begin{align}
fit2(x)=&5.49888\times10^{-9}(-1.67134\times10^6+93810.2x+1.73499\times10^6x^2-845361 x^3+\notag\\&1.3559\times10^6x^4+734284 x^5-184196 x^6-715654 x^7-502434 x^8)\times\notag\\&[-\tan^{-1}\big(1000000(x-1.00037)\big)+\tan^{-1}\big(1000000(x-0.999345)\big)]\times\notag\\&(1+\sin(1.76694\times10^7-3.53555\times10^7x+1.76861\times10^7x^2))\notag
\end{align}\]
另：产生切比雪夫多项式的命令：

1. Table[Cos[i ArcCos[x]] // TrigExpand, {i, 10}] // MatrixForm
   
2. Table[Sin[i ArcSin[x]] // TrigExpand, {i, 10}] // MatrixForm

复制代码

Jack315

elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨，而是刻意不动脑子不学习不动手的笨．

@Jack315 :  建议用切比雪夫理论一举得到支集 ...

使用切比雪夫多项式进行拟合的结果：


最大误差（大约）为：\(9.28362\times10^{-9}\)

拟合公式 (第一类 Chebyshev 多项式，n=50)：
\(a_0+\sum_{i=1}^{n}a_iChebyshevT(2i,x)\)

代码：
答案.rar (5.21 KB, 下载次数: 4)

2023-9-25 15:23 上传

点击文件名下载附件