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

笨笨

Jack315 发表于 2023-9-5 13:48
一个函数展开成幂级数时为：
\[\sum_{n=0}^{\infty}c_n(x-x_0)^n\]
在构造补偿函数 \(\varphi(x)\) 时，\(x ...

可以试试这个函数族,不知道怎么样:

Jack315

好，这个可以试试。

90#的这个式子 \(\sum_{n=0}^{\infty}c_n(x-x_0)^n\) 意味着在 \(x=x_0\) 处的误差为 \(c_0\) 。
而在 \(x=x_0\) 处的 \(c_0\) 就是被展开的误差函数值，即 \(c_0=err(x_0)=2F1(x_0)-Fit(x_0)\)。
所以绕了一圈回到了拟合问题经典方法上了。
一般的拟合问题是给出一组数据，即观察值，然后拟合一个函数。
这里不同的地方是观察值是由函数 \(2F1(x)\) 直接生成的，而不是通过实验得到的。

玩了这么久，小结 2 点：
1. 由于 \(2F1(x)\) 没有等价的初等函数表示，意味着由任何初等函数组合构成的拟合函数
都必然存在误差，也就是之前说的拟合精度有个“天花板”。
（理论上）只有分段的拟合精度才没有“天花板”。
2. 构造拟合函数应该从形状相似的函数开始，
通过修改拟合函数的形式并优化相应参数来提高拟合精度。

还是给个例子。
假设拟合函数的形式为：
\[f(x)=a+\frac{bx^2}{c+\sqrt{d+ex^2}}\]
用经典的拟合方法得到：
a=1.00001566169715
b=2.99565328063875
c=9.99609869674539
d=3.98258098897837
e=-3.03063398716602
此时的最大误差为：0.000191225 。
这个最大误差就是这个拟合函数形式精度的“天花板"。
（参考）拉马努金公式的最大误差为：0.000512272 。

LZ 搜寻补偿函数的方法实际上就是寻找
超几何函数与拉马努金公式之差的拟合函数。
画图出来看看：

1. LogLogPlot[
   
2. Hypergeometric2F1[-(1/2), -(1/2), 1, x^2] - (1 + (3 x^2)/(10 + Sqrt[4 - 3 x^2])), {x, 0, 1},
   
3. WorkingPrecision -> 100,
   
4. GridLines -> Automatic,
   
5. Frame -> True,
   
6. PlotRange -> All]

复制代码

图形几乎是一条直线，但在接近 x=1 的附近处开始往上翘起。
这提示我们函数形式或可以从 \(\ln(\varphi(x))=k\ln(x)+b\) 开始（修正），即：
\[\varphi(x)=e^{k\ln(x)+b}=e^bx^k\]

Jack315

笨笨 发表于 2023-8-18 20:57

18#,20# 的误差曲线与这个函数的形状类似，
看看能不能在这个函数上下点功夫：

1. Plot[0.000000182658599579 Sin[15.7 x^10], {x, 0, 1}, GridLines -> Automatic, Frame -> True]

复制代码

误差曲线上有 6 个零点，3 个极大值，2 个极小值。
用这 11 个点拟合 Sin 函数里的 10 次多项式（11 个参数），
不知道效果会如何。

Jack315

Jack315 发表于 2023-9-7 08:32
18#,20# 的误差曲线与这个函数的形状类似，
看看能不能在这个函数上下点功夫：

要把所有的峰都对齐，看样子又有点难度。
这里再给一个脱胎于 Weill 分布函数的一种函数形式：
\[c\bigg(\frac{1-bx}{a}\bigg)e^{-\big(\frac{1-bx}{a}\big)^2}\]

1. weibull[x_, {a_, b_, c_}] := c ((1 - b x)/a) E^-((1 - b x)/a)^2
   
2. 
   
3. Plot[{weibull[x, {0.2, 2, 1}], weibull[x, {0.2, 1, 1}]}, {x, 0, 1},
   
4. PlotRange -> All,
   
5. GridLines -> Automatic,
   
6. Frame -> True]

复制代码

一个这样的函数对齐两个峰。也可以对齐最后一个峰。
这样或许可以将精度提高到 1e-8 以下。

Jack315

Jack315 发表于 2023-9-7 17:22
要把所有的峰都对齐，看样子又有点难度。
这里再给一个脱胎于 Weill 分布函数的一种函数形式：
\[c\bigg( ...

这个函数形式能行，这是对前两个峰补偿的效果：

一轮降峰应该能到 1e-8 以下。
多轮降峰最后会不会也有“天花板”还要试了才知道。

Jack315

对于剩余误差 \(F(x)-G(x)-H(x)\) 再进行补偿的设想。

剩余误差有五个极值点，六个零点。
补偿的关键在于五个极值点的位置 (x) 必须精确对准。
通常的基本初等函数组合，如幂函数的和，很难做到这一点。
为避免使用分段函数，使用由正态分布函数导出的一个补偿函数：
\[\frac{a}{s}e^{-b(\frac{x-m}{s})^2}\]
让五个这样的补偿函数分别对准五个极值点。

补偿的方法由两种：
1. 由五个补偿函数的和拟合剩余误差，即最小化下列积分：
\[\int_0^1\bigg[err(x)-\sum_{i=1}^5\frac{a_i}{s_i}e^{-b_i(\frac{x-m_i}{s_i})^2}\bigg]^2dx\]
2. 由五个补偿函数的和在极值点精确拟合角度，而在即它点上曲线尽量保持顺滑。
即最小化下列和式（应该为零，可能还需要其它优化约束条件）：
\[\sum_{i=1}^5\bigg[\theta_i-\frac{a_i}{s_i}e^{-b_i(\frac{x_i-m_i}{s_i})^2}\bigg]^2\]
其中：\(x_i\) 为五个极值点位置，\(\theta_i=n\frac{\pi}{2}\)，为极值对应的角度。
然后用 \(k\sin[\sum_{i=1}^5\frac{a_i}{s_i}e^{-b_i(\frac{x_i-m_i}{s_i})^2}]\) 函数来拟合误差函数 \(err(x)\) 。

第一种方法比较有望能实现，精度应该也能到 1e-8 以下；
第二种方法若能实现，精度有望达到更高的水平。

笨笨

这个函数形式能行，这是对前两个峰补偿的效果：

请问这个函数形式是啥?

笨笨

Mathematica具有“滑动条”功能

1. Manipulate[
   
2. Plot[{(Hypergeometric2F1[-(1/2), -(1/2), 1, x^2] - 1 - (3 x^2)/(
   
3.      10 + Sqrt[4 - 3 x^2]) -
   
4.      3/2^17 x^10 ((4/\[Pi] - 14/11) 2^17/3)^(
   
5.       x^2/(1 +
   
6.         425.524154 x^-2.4565895 (1 -
   
7.            x^0.04701004)^0.94801807)^0.1537177)),
   
8.    a Sin[b x^c + d x + e], -1.3830579792539766`*^-7,
   
9.    1.3830579792539766`*^-7}, {x, 0, 1}], {a, -20, 20}, {b, -20,
   
10.   20}, {c, -20, 20}, {d, -20, 20}, {e, -20, 20}]

复制代码

elim

还没适应本论坛的 MathJax 设置。所以贴个图：

Jack315

elim 发表于 2023-9-9 23:28
还没适应本论坛的 MathJax 设置。所以贴个图：

差不多想到一个方向上了……

设拟合函数为：\(f_{Fit}(x;a,b,c)=a+bx+cx^2\)；
窗函数为：\(f_{Win}(x)=e^{-(\frac{x-\mu}{\sigma})^2}\)
则两函数的积 \(f(x;a,b,c)=f_{Fit}(x;a,b,c) \times f_{Win}(x)\)
就可以得到某个区间内的拟合函数。而区间外接近为零。
函数集：\(\sum_{i=1}^nf(x;a_i,b_i,c_i)\) 即可将任意函数拟合到任意想要的精度。
其实这个就是分段拟合的思路，则是不出现 \(Picewise(x;a,b)\) 函数而已。

如果拟合函数使用更多项：\(f_{Fit}(x;\lambda_0,\lambda_1,...\lambda_k)=\sum_{i=0}^k\lambda_ix^i\);
而窗函数使用：\(f_{Win}(x)=\bigg\lfloor\frac{|1+x|}{1+|x|}\bigg\rfloor=\begin{cases}1&\text{if }x\ge0\\0&\text{if }x<0\end{cases}\)
则拟合效果更好，逼近速度也更快。