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

笨笨

Jack315 发表于 2023-9-2 10:29
【分享拉马努金 (Ramanujan) 公式推导过程】

首先将超几何函数展开：

笨笨

Jack315 发表于 2023-9-2 10:47
【补偿函数的构造】
画出误差函数的曲线：

所拟合函数的定义域是0<x<1，而正割函数在\(\left( 0，\frac{\pi}{2}\right)\)图像确实相仿，但它们定义域不同，还要满足我楼上刚发图片上的要求，想想怎么转化之间关系。

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\) 处依然存在，但不再是实数？
还是说软件在照章办事？