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

Jack315

为了表示对大师的敬意，也为了尊重 LZ 出的题，
还是希望能得到一个答案。先来个临时结果：
函数拟合误差优化.rar (54.3 KB, 下载次数: 1)

2023-8-26 10:55 上传

点击文件名下载附件

对这个临时结果做个简单说明，细节部分看文档。
在 28# 提了个思路：

【解题思路】
假设给一组 (a,b,c,d) 的值，利用 (EQ:4) 定义的
Err 函数就能求出 x 在其定义域 [0,1] 范围内的最大值。
将这个过程包装到一个新的函数中：

这个新的函数在文档中被称为模型的”评价函数“。
在数学（优化）的术语中，称之为响应函数。

在文档中虽然误差函数都一样，但使用了两个不同的响应函数。
一个是求误差函数的绝对值的最大值；
另一个是求误差函数平方的积分。由于积分运算比较费时，
这个响应函数又改成了离散化后的误差平方和的平均值。
两者精度接近，但计算速度提高很多。

第一个响应函数运行了很长时间，CPU 也开始冒烟了，
但仍然还没有结果出来。原因应该是响应曲面不光滑，
可能会是比较粗糙：参看 20# 的图。
NMinimize 命令是求全局最优解的，因而导致了计算量剧增。

第二个响应函数瞬间就出结果了：(a->0, b->0, c->0, d->0) 。
有一定的补偿效果，看起来至少不是不靠谱了。
由于在 NMinimize 命令的设置里，只使用了参数大于零的空间 (1/16) 。
而搜索结果指向这部分空间的边界，说明这部分空间里没有极值，
或者说最优解至少应该有一个参数必须小于零。

补偿函数的计算只是在 x=0 这个点上有问题，后续代码再作修改。

笨笨

Jack315 发表于 2023-8-26 11:25
为了表示对大师的敬意，也为了尊重 LZ 出的题，
还是希望能得到一个答案。先来个临时结果：

本楼作为补充说明用：

昨晚整理的，这是四参数拟合的效果图看看有没有可靠的信息。



另外补充一下，一个小细节转换：

Jack315

能分享下这个公式的由来吗？
\[\begin{matrix}H(x,a,b,c,d)=\frac{3x^{10}}{2^{17}}(1+\frac{[\frac{79}{48}+(\mu-1-\frac{79}{48})x]x^2}{[1+x^a(1-x^b)^c]^d}),&\mu=\frac{2^{17}}{3}(\frac{4}{\pi}-\frac{14}{11})\end{matrix}\tag{EQ:3}\]

Jack315

这个题有可能会走进死胡同了……

先用下列代码定义相关的函数：

1. ClearAll["Global`*"]
   
2. x0 = x -> 0.997;
   
3. param = {a -> 0.53, b -> 2.75, c -> 0.82, d -> 3.38};
   
4. target[x_] := Hypergeometric2F1[-(1/2), -(1/2), 1, x^2]
   
5. ramanujan[x_] := 1 + (3 x^2)/(10 + Sqrt[4 - 3 x^2])
   
6. compensation[x_, a_, b_, c_, d_] := ((3 x^10)/2^17
   
7.   (1 + ((79/48 + (\[Mu] - 1 - 79/48) x) x^2) /
   
8.   (1 + x^a (1 - x^b)^c)^d)) /. \[Mu] -> 2^17/3 (4/\[Pi] - 14/11)
   
9. err[x_, a_, b_, c_, d_] := target[x] - ramanujan[x] - compensation[x, a, b, c, d]

复制代码

然后逐一画出与各参数对应的误差函数的图来：

1. Plot3D[{err[x, varA, b, c, d] /. param, 0},
   
2. {x, 0, 1}, {varA, 0, 2 a /. param},
   
3. WorkingPrecision -> 30,
   
4. AxesLabel -> Automatic,
   
5. PlotLegends -> {"误差", "零平面"}]

复制代码

1. Plot3D[{err[x, a, varB, c, d] /. param, 0},
   
2. {x, 0, 1}, {varB, 0, 2 b /. param},
   
3. WorkingPrecision -> 30,
   
4. AxesLabel -> Automatic,
   
5. PlotLegends -> {"误差", "零平面"}]

复制代码

1. Plot3D[{err[x, a, b, varC, d] /. param, 0},
   
2. {x, 0, 1}, {varC, 0, 2 c /. param},
   
3. WorkingPrecision -> 30,
   
4. AxesLabel -> Automatic,
   
5. PlotLegends -> {"误差", "零平面"}]

复制代码

1. Plot3D[{err[x, a, b, c, varD] /. param, 0},
   
2. {x, 0, 1}, {varD, 0, 2 d /. param},
   
3. WorkingPrecision -> 30,
   
4. AxesLabel -> Automatic,
   
5. PlotLegends -> {"误差", "零平面"}]

复制代码

这些图说明：
1. 误差函数的过零点（解）就在 param 这个点附近；
2. 在过零点附近，误差函数的曲面如陡崖峭壁——斜率非常大。

这意味着四个参数要精确定位到零点，其精度要求非常高。
是非常非常非常的高。而导致这种情况的原因
很可能就是补偿函数的形式导致了误差函数不是单调的。
或许应该重新寻找合适的拟合函数形式。
包括大师提供的这个拟合函数。
当然也可以继续保留——下一个贴说明原因。

Jack315

【寻找拟合函数的迭代过程】
任意给定一函数 \(\begin{matrix} y=F_0(x),&x\in[a,b]\end{matrix}\)，
将这个函数看作是与目标函数 \(y=0\) 的误差函数。
目标是寻找一组基本初等函数 \(\begin{matrix}y_i=F_i(x),&i=1,2,...,n\end{matrix}\) 的组合，
以使叠加上这组函数后，误差函数与横坐标之间的面积最小。

将这些基本初等函数加上横、竖两个方向进行伸缩和移位的参数，
并与这些函数某些特定的组合函数作为备选函数，按下列步骤进行迭代：
1. 从备选函数中选出与误差函数形状最接近的一个函数。
2. 优化这个函数的参数。
3. 叠加上这个函数后形成新的误差函数。
4. 判断拟合精度是否满足要求，
不满足要求就重复步骤 1~3 的迭代过程。

这是一个设想，直觉上应该可行。
因为基本初等函数大部分都有比较好的特性，
只要每次叠加上选出的函数后，
误差函数与坐标轴之间的面积都能（大幅）减少就行。

所以实际过程是将每个备选函数的参数先进行优化，
与误差函数之间所围面积最小的就是被选中的函数。

在选择备选函数及其相应参数的选择上，
应避免导致计算困难的形式，以便进行参数的优化。
大师提供的拟合函数也符合这些要求。

笨笨

Jack315 发表于 2023-8-27 08:41
【寻找拟合函数的迭代过程】
任意给定一函数 \(\begin{matrix} y=F_0(x),&x\in[a,b]\end{matrix}\)，
将 ...

用简单的初等函数构造拟合函数貌似行不通？

不过楼上的elim老师曾经说过：

现在可设法用“初等超越函数”\(\varphi\)拟合楼上的 \(Er\), 使 \(\|E_r-\varphi\|_{\infty}\le 10^{-12}\)

作为阶段性的有理论合实际意义的结果.

附：论坛处于维护中，还是无法上传图片。

Jack315

【拟合函数相关定义】

备选初等函数清单，不包括其组合：
1) 常数函数：\(f(x)=c\)
2) 幂函数：\(f(x)=x^a\)
3) 指数函数：\(f(x)=a^x\)
4) 对数函数：\(f(x)=\log_a(x)\)
5) 三角函数：\(f(x)=\sin(x),\ \cos(x),\ \tan(x),\ \cot(x),\ \sec(x),\ \csc(x)\)
6) 反三角函数：\(f(x)=\sin^{-1}(x),\ \cos^{-1}(x),\ \tan^{-1}(x),\ \cot^{-1}(x),\ \sec^{-1}(x),\ \csc^{-1}(x)\)

进行组合的运算：
（有限次的）加、减、乘、除、乘方、开方。

笨笨

Jack315 发表于 2023-8-26 16:50
能分享下这个公式的由来吗？
\[\begin{matrix}H(x,a,b,c,d)=\frac{3x^{10}}{2^{17}}(1+\frac{[\frac{79}{48 ...

原构造思想来源于楼上“elim”老师的贴子，楼主只是作小幅度改动。不过大致构造思想如下：



另外楼主这几天重新构造一个，效果不错。

\(Er\left( x \right) ≈ {}_2{F_1}\left( { - \frac{1}{2}, - \frac{1}{2},1,{x^2}} \right) - 1 - \frac{{3{x^2}}}{{10 + \sqrt {4 - 3{x^2}} }} - \frac{3}{{{2^{17}}}}{x^{10}}{\mu ^{\frac{{{x^2}}}{{{{\left( {1 + 425.524154{x^{ - 2.4565895}}{{\left( {1 - {x^{0.04701004}}} \right)}^{0.94801807}}} \right)}^{0.1537177}}}}}}\)

误差精度：\(1.383057979\cdots\times10^{-7}\)

Jack315

笨笨 发表于 2023-8-28 11:02
原构造思想来源于楼上“elim”老师的贴子，楼主只是作小幅度改动。不过大致构造思想如下：

感谢 LZ 分享！

这是补偿函数一步到位的思路。
无论最终用哪一类函数族作为 \(\varphi(x)\)，都可能会导致下列问题：
如果误差函数比较粗糙，就是坑坑洼洼的那种，
如 \(f(x,y)=\sin(x^2+y^2)\) 这样的曲面，
软件命令就可能直接无法找到全局最优解，除非运气好。
如果是用光滑函数进行迭代，可能软件命令就能直接搞定。

为了解决找到可能的全局最优解这个问题，或可考虑先进行一个初步的搜索。
期望这个搜索能找到至少在平均意义上拟合效果最佳的区域。
这个就是 elim 在 21# 提出的思想。
现在正在尝试实现这个想法，对一步到位的这类函数族作最后一次努力。

在找到这个最佳区域后，再用软件的命令在这个区域内找出局部最优解。
希望这次能有令人满意的结果。

笨笨

Jack315 发表于 2023-8-28 12:21
感谢 LZ 分享！

这是补偿函数一步到位的思路。

33楼elim的方法有望突破10^(-12).



设拉马努金拟合误差函数为\(\psi_0,\;\psi_n\)是函数\(Er_{n-1}:=\psi_0-\cdots-\psi_{n-1}\)
的拟合， 使得 \(\lVert Er_k\rVert /\lVert Er_{k-1}\rVert\le 10^{-2},\:\lVert f\rVert
:=\max|f|((0,1))\)
则对某\(n\le 6\) 就有 \(\lVert Er\rVert \le 10^{-12}\).

以上方法可称为辗转拟合法．



相信这里的探讨的理论阐述会相当繁琐，或许可以从优化参数序列的(包络)极限来理解？总之很妙很有趣！

现在可设法用初等超越函数\(\varphi_2\)拟合楼上的 \(Er\), 使 \(\|E_r-\varphi_2\|_{\infty}\le 10^{-12}\)
作为阶段性的有理论合实际意义的结果.