找回密码
 欢迎注册
楼主: 笨笨

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

  [复制链接]
 楼主| 发表于 2023-9-14 22:07:21 | 显示全部楼层
本帖最后由 笨笨 于 2023-9-14 22:46 编辑

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


\(\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 \)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-15 15:09:51 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-15 17:55 编辑
笨笨 发表于 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. 误差函数为经部分拟合后的剩余误差。

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

各位路过的高手多提宝贵意见,谢过先
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-9-16 18:31:22 | 显示全部楼层
本帖最后由 笨笨 于 2023-9-16 18:36 编辑
Jack315 发表于 2023-9-15 15:09
敢问 LZ 是什么专业 。言归正传……

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


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

86楼是楼主的个人极限,希望志同道合者有人能突破这个极限。最好对拉马努金公式再拟合找到一个简洁的再拟合公式一举突破10^-12,那就牛逼了!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-20 16:01:17 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-20 16:13 编辑
笨笨 发表于 2023-8-30 17:20
这是楼主目前找的比较不错的参数与拟合函数族,小结一下:

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


这是在 18# 前辈的基础上用 Sin 函数一次性再补偿的效果:
再补偿误差_1.png
再补偿误差_2.png
再补偿误差_3.png
再补偿误差_4.png

补偿函数为:
\[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\) 的区域。
最终仍然需要用加窗函数的办法来解决。

单峰函数和窗函数有比较多的细节问题……有时间再接着玩。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-9-20 20:14:59 | 显示全部楼层
Jack315 发表于 2023-9-20 16:01
这是在 18# 前辈的基础上用 Sin 函数一次性再补偿的效果:

看看86楼我的历史记录,能否突破?

点评

期待,感谢先生  发表于 2023-9-20 21:50
等加窗的情况搞清楚后,突破最佳纪录应该有希望的。  发表于 2023-9-20 21:35
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-9-21 22:30:55 | 显示全部楼层
elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨,而是刻意不动脑子不学习不动手的笨.

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

请问怎么用切比雪夫理论一举得到最佳多项式一致逼近?可否具体点????
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-22 02:37:01 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-22 04:23 编辑
笨笨 发表于 2023-9-20 20:14
看看86楼我的历史记录,能否突破?


稍微再改进了一点点:
再补偿误差_5.png

最终误差公式:
\[\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
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-25 15:23:34 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-25 16:42 编辑
elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨,而是刻意不动脑子不学习不动手的笨.

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


使用切比雪夫多项式进行拟合的结果:
拟合误差.png
最大误差.png
最大误差(大约)为:\(9.28362\times10^{-9}\)

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

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

点评

两位前辈功底深厚,点个赞!Ickiverar 的代码优雅、高效,值得潜心学习。  发表于 2023-9-25 15:33
这个题玩得比较久,学到了不少知识,包括数学和编程。也荒废了不少正事 :(  发表于 2023-9-25 15:30
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-29 00:16:35 | 显示全部楼层
我一直想用级数展开继续搞这个,但是也不好搞出形式简洁而且误差较低的。

无限增加级数项,倒是可以无限降低误差,但形式会以项数的平方变复杂,非常糟糕。

ErApprox.zip (55.81 KB, 下载次数: 3)

评分

参与人数 1威望 +12 金币 +12 贡献 +12 经验 +12 鲜花 +12 收起 理由
uk702 + 12 + 12 + 12 + 12 + 12

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-29 11:24:18 | 显示全部楼层
Ickiverar 发表于 2023-9-29 00:16
我一直想用级数展开继续搞这个,但是也不好搞出形式简洁而且误差较低的。

无限增加级数项,倒是可以无限降 ...

这个公式的结果肯定不是最好,个人觉得相对比较简单,可能也相对容易扩展。
2023-09-29_112057.png
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-10-14 22:12 , Processed in 0.029015 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表