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

Jack315

笨笨 发表于 2023-9-2 20:45

除了拉马努金公式，应该还有其它的函数形式可以达到类似的效果，
只要是包含4个参数，就可以抵消级数展开式中的4个项……
题目中拉马努金公式是必须的吗？

【设想】
除了用一个函数进行补偿外，或许还可能先将拟合函数旋转一个角度，
然后再进行伸缩，使得拟合函数与2F1函数在定义域两端完全重合。
这样误差就可能会减少一（大）部分。这时再用一个函数去补偿或许更容易一些。

笨笨

Jack315 发表于 2023-9-2 21:56
除了拉马努金公式，应该还有其它的函数形式可以达到类似的效果，
只要是包含4个参数，就可以抵消级数展开 ...

按照你的思路走，可否构造一个函数族来且得到什么样的公式？

Jack315

笨笨 发表于 2023-9-2 22:29
按照你的思路走，可否构造一个函数族来且得到什么样的公式？

按照设想先选出单纯由初等基本函数构成的拟合函数：
\[\begin{align*}&1+\frac{3x2}{10+\sqrt{4-3x^2}}\tag{1}\\&1+\frac{1}{4}x^2\tag{2}\\&\frac{64-3x^4}{64-16x^2}\tag{3}\\&-1+2e^{\frac{1}{8}x^2}\tag{4}\\&\frac{1}{3}(1+2\cosh(\frac{\sqrt{3}}{2}x))\tag{5}\\&1+\frac{8}{3}\tan^2(\frac{1}{4}\sqrt{\frac{3}{2}}x)\tag{6}\\&\frac{1}{3}\bigg[-7+10\sec\bigg(\frac{1}{2}\sqrt{\frac{3}{5}}x\bigg)\bigg]\tag{7}\\&1+\frac{4}{3}\arcsin^2(\frac{\sqrt{3}}{4}x)\tag{8}\end{align*}\]
(1) 式为拉马努金公式作为参考。
到目前为止构造出的拟合函数，都是单调上升的，
因此最大误差都发生在 \(x=1\) 处。
比较结果拉马努金公式的误差最小，形式也比较优雅。

经过旋转拉伸后，公式表达形式可能会比较复杂。
而由三角函数构成的拟合函数在这一点上或许有一定的优势。
反正 (1) 式经这样处理后，出来的公式特别复杂。

接下来先看看误差能减少多少。
如果再需要补偿的时候，误差函数都是两头为零、中间拱起的形状。
这种情形找补偿函数应该会容易一些。

附件是构造拟合函数的代码，供参考：
超几何函数拟合函数.rar (27.87 KB, 下载次数: 2)

2023-9-3 10:01 上传

点击文件名下载附件

各拟合函数的性能：


拟合函数的上升速度都比 2F1 函数慢。
如果有一个超过的话，一个线性组全就能使两头端点重合了。

笨笨

Jack315 发表于 2023-9-3 10:03
按照设想先选出单纯由初等基本函数构成的拟合函数：
\[\begin{align*}&1+\frac{3x2}{10+\sqrt{4-3x^2}}\t ...

可以看出拉马努金的公式：\(1+\frac{3 x^2}{10+\sqrt{4-3 x^2}}\)与椭圆周长原级数公式重合度较高，因此只需对拉马努金公式再拟合即可。

所以现在的目标就是找到更好的\(\varphi \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}} }} + \varphi \left( x \right)\)

下面是楼主找到的:
\(Er\left( x \right) \approx {}_2{F_1}\left( { - \frac{1}{2}, - \frac{1}{2},1,{x^2}} \right) - \left( {1 + \frac{{3{x^2}}}{{10 + \sqrt {4 - 3{x^2}} }}} \right) - \frac{3}{{{2^{17}}}}{x^{10}}{\left[ {\left( {\frac{4}{\pi } - \frac{{14}}{{11}}} \right)\frac{{{2^{17}}}}{3}} \right]^{\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}\)

你能找到更好的\(\varphi \left( x \right)\)吗?

Jack315

问题主要在 x=1 这个奇点及其邻域。
虽然不是剧烈的振荡，但函数的各阶导数却是快速上升。
将这个 2F1 函数的各阶导数画出来就可以看出
要用简单地初等函数获得高精度的拟合有多困难。
下面给个例子。

拟合区间为[0.999, 1.000]：

1. rng = {0.999, 1.000}

复制代码

拟合函数为二次函数：

1. f[x_, {a_, b_, c_}] := a + b x + c x^2

复制代码

用于 Table 命令中的计算 x 值的辅助转换函数：

1. t[n_, {xmin_, xmax_}] := n/2 (xmax - xmin) + xmin

复制代码

利用三个点上两函数值相等来求出 (a,b,c) ：

1. sol = N[Solve[Table[f[t[i, rng], {a, b, c}] ==
   
2.      Hypergeometric2F1[-(1/2), -(1/2), 1, (t[i, rng])^2], {i, 0, 2}],
   
3.    {a, b, c}], 100]

复制代码

最后画出误差曲线。

1. Plot[Hypergeometric2F1[-(1/2), -(1/2), 1, x^2] - f[x, {a, b, c}] /. sol, {x, rng[[1]], rng[[2]]},
   
2. WorkingPrecision -> 100,
   
3. GridLines -> Automatic,
   
4. Frame -> True]

复制代码

可以看出精度大致在 1E-8 到 1E-9。

所以理论上，只要拟合区间足够小，
就能达到任意想要的精度，而且是【有限】项的和。
也就是分段拟合能获得任意高的精度。
在实际上，则有可能会受计算精度所限。

而回过头来看楼主找的这个补偿函数。
能达到 1E-7 的水平已经非常优秀。
因为从拟合的角度而言，一个拟合函数要比分段的多个拟合函数好，
只要能达到所需的精度要求就行。

笨笨

Jack315 发表于 2023-9-5 08:26
问题主要在 x=1 这个奇点及其邻域。
虽然不是剧烈的振荡，但函数的各阶导数却是快速上升。
将这个 2F1 函数 ...

能否无需分段且对拉马努金公式再拟合达到10^（-8）或10^（-12）？
意思就是说能不能找到比楼主更好的拉马努金公式再拟合函数族？

Jack315

笨笨 发表于 2023-9-5 10:21
能否无需分段且对拉马努金公式再拟合达到10^（-8）或10^（-12）？
意思就是说能不能找到比楼主更好的拉马 ...

LZ 构造的补偿函数应该继续发扬光大……真想一睹 LZ 灵感的光亮。
以这个思路下去到 1E-8 以下的精度非常有可能。
精度要求再高，理论上不是没有可能，
但恐怕会遇到计算上的问题，包括：
1. 计算精度
2. 计算时长
让人想起来人家那种挖矿，挖呀挖呀挖……LZ 这个也能挖到矿吗？

Jack315

一个函数展开成幂级数时为：
\[\sum_{n=0}^{\infty}c_n(x-x_0)^n\]
在构造补偿函数 \(\varphi(x)\) 时，\(x_0\) 应尽可能选择接近 1 。
这样构造出来的补偿函数，能使 \(x=x_0\) 处的拟合误差最小
（\(c_0=2F1(x_0)\) 时为零），而 \(x=x_0\) 的领域内，误差也会比较小。
由于 \(x=1\) 处是奇点，\(x_0\) 无法选择 1 ，否则这是最佳选择。
基于这一点，LZ 的补偿函数有更新吗？

如果选择 \(x_0=1\) ，函数 2F1 的级数展开式中的系数出现复数。
不是太懂这是什么意思，有人能帮解释一下吗？
这是说函数 2F1 的各阶导数在  \(x=1\) 处依然存在，但不再是实数？
还是说软件在照章办事？

笨笨

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\]