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

elim

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

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

补充内容 (2023-9-15 05:35):
应为非初等函数。

笨笨

elim 发表于 2023-8-23 02:54
相信这里的探讨的理论阐述会相当繁琐，或许可以从优化参数序列的(包络)极限来理解？总之很妙很有趣！

现 ...

10^(-12)是楼主一直以来的目标，这样的初等函数族或初等超越函数族估计没人能找出来！

光搞个10^（-7）都这么费劲，10^（-12）想都别想！

Jack315

【被拟合函数】
有一个比较复杂的函数：
\[\begin{matrix} F(x)=\sum_{n=0}^ \infty \begin{pmatrix}1/2\\n\end{pmatrix}^2x^{2n},&x\in[0,1]\end{matrix}\]
设想找一个由初等函数构成的函数 \(G(x)\) 来拟合 \(F(x)\)，
使得两者之间的误差函数 \(Er(x)=F(x)-G(x)\) 的绝对值在定义域范围内越小越好。
目前希望能做到  \(|Er(x)|_{MAX}<10^{-8}\) 的水平 。

在题目中给出了这样一个等式：
\[\begin{pmatrix}1/2\\k\end{pmatrix}=\frac{1}{k!}\prod_{j=0}^{k-1}\bigg(\frac{1}{2}-j\bigg)\]
这个等式只有在 k 取正整数时成立。

现在主要关注等式的左边，因为它就出现在 \(F(x)\) 的定义中。
在 Mathematica 中它对应 \(Binomial[\frac{1}{2},k]\) 函数。
用下列命令简化 \(F(x)\) 函数：

1. FullSimplify[\!\(TraditionalForm\`
   
2. \*UnderoverscriptBox[\(\[Sum]\), \(n = 0\), \(\[Infinity]\)]\((
   
3. \*SuperscriptBox[\(Binomial[
   
4. \*FractionBox[\(1\), \(2\)], n]\), \(2\)]
   
5. \*SuperscriptBox[\(x\), \(2  n\)])\)\)]

复制代码

\[\text{FullSimplify}\left[\sum _{n=0}^{\infty } \binom{\frac{1}{2}}{n}^2 x^{2 n}\right]\]
结果得到 \(Hypergeometric2F1[-\frac{1}{2},-\frac{1}{2},1,x^2]\)。
所以定义的这个比较复杂的函数实际上是一个超几何函数：
\[F(x)=Hypergeometric2F1[-\frac{1}{2},-\frac{1}{2},1,x^2]\tag{EQ:1}\]

当然，如果用上述等式右边来定义 \(F(x)\)
可以得到另外一个等价的函数。这里不作进一步讨论。

在 Mathematica 中用下列命令定义 \(F(x)\) ：

1. f[x_] := Hypergeometric2F1[-(1/2), -(1/2), 1, x^2]

复制代码

Jack315

【拟合函数】
在题目中给出了这样一个拟合函数：
\[G(x)=1+\frac{3x^2}{10+\sqrt{4-3x^2}}\tag{EQ:2}\]
在 Mathematica 中用下列命令定义 \(G(x)\) 函数：

1. g[x_] := 1 + (3 x^2)/(10 + Sqrt[4 - 3 x^2])

复制代码

接下来用下列命令画出图来看一下这两个函数之间的误差曲线：

1. LogPlot[(f[x] - g[x]), {x, 0, 1},
   
2. PlotRange -> All,
   
3. GridLines -> Automatic,
   
4. Frame -> True,
   
5. PlotLabel -> "待拟合函数与拟合函数的误差"]

复制代码

说明：图能不上不上了，新人权限太低，一直被限制。
大家可以自己在软件里画出来看。

结果发现，随着 x 的增加，误差快速增大，使得拟合精度变差。
所以就需要再用一个修正函数 \(H(x)\) 来修正这个误差。

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

用下列命令定义 \(H(x)\) 函数：

1. h[x_, a_, b_, c_, d_] := ((3 x^10)/2^17 (1 + ((79/48 + (\[Mu] - 1 - 79/48) x) x^2)/(1 + x^a (1 - x^b)^c)^d)) /. \[Mu] -> 2^17/3 (4/\[Pi] - 14/11)

复制代码

题目要求能找出这样一组 \((a,b,c,d)\) 的值，
使得拟合的糖度能达到要求。

Jack315

【误差函数】
定义下列拟合的误差函数：
\[Err(x,a,b,c,d)=F(x)-G(x)-H(x,a,b,c,d)\tag{EQ:4}\]

用下列命令定义误差函数：

1. err[x_, a_, b_, c_, d_] := f[x] - g[x] - h[x, a, b, c, d]

复制代码

题目中给了个例子：
例：当 \((x,a,b,c,d)=(0.997,0.53,2.75,0.82,3.38)\) 时，\(Er\approx6\times10^{-6}\)。

可以使用下列命令来查看误差函数的情况：

1. (*定义参数*)
   
2. x0 = x -> 0.997;
   
3. param = {a -> 0.53, b -> 2.75, c -> 0.82, d -> 3.38};
   
4. (*查看误差函数值*)
   
5. err[x, a, b, c, d] /. x0 /. param
   
6. (*查看定义域内误差曲线*)
   
7. Table[Plot[err[x, a, b, c, d] /. param, {x, (i - 1)/4, i/4},
   
8.   PlotRange -> All,
   
9.   GridLines -> Automatic,
   
10.   PlotLegends -> "Expressions",
    
11.   Frame -> True,
    
12. {i, 4}]

复制代码

笨笨

Jack315 发表于 2023-8-24 17:02
【误差函数】
定义下列拟合的误差函数：
\[Err(x,a,b,c,d)=F(x)-G(x)-H(x,a,b,c,d)\tag{EQ:4}\]

感谢先生闲暇时间回帖.有以下内容可供参考。

Jack315

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

1. maxErr[aVal_, bVal_, cVal_, dVal_] :=
   
2. Module[{a = aVal, b = bVal, c = cVal, d = dVal},
   
3.   MaxValue[{Abs[err[a, b, c, d, x]], 0 <= x <= 1}, x \[Element] Reals]
   
4.   ]

复制代码

然后用下列命令进行求解：

1. sol = NMinimize[{maxErr[a, b, c, d], {a, b, c, d} \[Element] Reals}, {a, b, c, d}]

复制代码

运行一段时间后，给出了结果……

结果不贴出来了，因为答案是错误的。
原因是 \(G(x)+H(x)\) 函数不是在整个参数空间都能计算出结果，
因而导致（自动的）优化过程步入了歧途。

解决的方法大致有两种：
1. 修改代码，使得代码能规避上述导致错误的原因……期待高人出手。
2. 利用相关的方法（如 DOE、RSM：即梯度优化），一步步的进行手工优化。
这个方法效率特别低，所以不建议，但最终应该能得到一个结果。
至少可以得到一个局部解。

笨笨

Jack315 发表于 2023-8-24 17:02
【误差函数】
定义下列拟合的误差函数：
\[Err(x,a,b,c,d)=F(x)-G(x)-H(x,a,b,c,d)\tag{EQ:4}\]

先生你好，代码少一个中括号：

Jack315

笨笨 发表于 2023-8-24 18:05
先生你好，代码少一个中括号：

复制的时候漏掉了。重新来过：

1. (*定义参数*)
   
2. x0 = x -> 0.997;
   
3. param = {a -> 0.53, b -> 2.75, c -> 0.82, d -> 3.38};
   
4. (*查看误差函数值*)
   
5. err[x, a, b, c, d] /. x0 /. param
   
6. (*查看定义域内误差曲线*)
   
7. Table[Plot[err[x, a, b, c, d] /. param, {x, (i - 1)/4, i/4},
   
8. PlotRange -> All,
   
9. GridLines -> Automatic,
   
10. PlotLegends -> "Expressions",
    
11. Frame -> True,
    
12. PlotLabel -> "超几何函数与（修正后）拟合函数的误差"],
    
13. {i, 4}]

复制代码