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

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}\)
则拟合效果更好，逼近速度也更快。

笨笨

Jack315 发表于 2023-9-10 02:55
差不多想到一个方向上了……

设拟合函数为：\(f_{Fit}(x;a,b,c)=a+bx+cx^2\)；

要考虑到定义域：0≤x≤1

Jack315

几种窗函数：
\[\begin{align*}&e^{-(\frac{x-m}{s})^2}\tag{1}\\&\frac{e^{2kx}(e^{2ak}+e^{2bk})}{(e^{2ak}+e^{2kx})(e^{2bk}+e^{2kx})}\tag{2}\\&\frac{e^{2kx}(-e^{2ak}+e^{2bk})}{(e^{2ak}+e^{2kx})(e^{2bk}+e^{2kx})}\tag{3}\\&\frac{1}{2tan^{-1}(\frac{b-a}{2}k)}\big(tan^{-1}[k(x-a)]-tan^{-1}[k(x-b)]\big)\tag{4}\end{align*}\]

(1) 式由正态分布概率密度函数导出；
(2)、(3) 式由双曲正切函数导出，
前者是两个双曲正切函数的乘积，
后者是两个双曲正切函数的差，
(4) 式是两个反正切函数的差。
顺便说一下：
1. 由 Weibull 分布概率密度函数导出的形式，感觉不太合适，暂时不考虑。
2. 两个双曲正切还可以构造出倾斜峰的形状，或可直接用于进行拟合。
3. 两个反正切函数的乘积的曲线形状与误差函数 \(2F1(x)-(1+\frac{3x^2}{10+\sqrt{4-3x^2}})\) 相似。

下列代码显示倾斜峰的曲线：

1. Plot[(Tanh[10 x - 5] Tanh[-20 x + 6] + 1), {x, 0, 1},
   
2. PlotRange -> All,
   
3. Filling -> Axis,
   
4. GridLines -> Automatic,
   
5. Frame -> True,
   
6. PlotLabel -> "倾斜峰函数"]

复制代码

下列代码定义这四个窗函数：

1. fwin0[x_, {m_, s_}] := E^-((x - m)/s)^2
   
2. fwin1[x_, {k_, a_, b_}] := (E^(2 k x) (E^(2 a k) + E^(2 b k)))/((E^(2 a k) + E^(2 k x)) (E^(2 b k) + E^(2 k x)))
   
3. fwin2[x_, {k_, a_, b_}] := (E^(2 k x) (-E^(2 a k) + E^(2 b k)))/((E^(2 a k) + E^(2 k x)) (E^(2 b k) + E^(2 k x)))
   
4. fwin3[x_, {k_, a_, b_}] := 1/(2 ArcTan[(b - a)/2 k]) (ArcTan[k (x - a)] - ArcTan[k (x - b)])

复制代码

下列代码显示窗函数曲线：

1. Plot[{fwin0[x, {0.5, 0.05}], fwin1[x, {50, 0.1, 0.3}],
   
2.   fwin2[x, {100, 0.4, 0.6}], fwin3[x, {1000, 0.7, 0.9}]}, {x, 0, 1},
   
3. PlotRange -> All,
   
4. Filling -> Axis,
   
5. GridLines -> Automatic,
   
6. Frame -> True,
   
7. PlotLabel -> "几种窗函数",
   
8. PlotLegends -> "Expressions"]

复制代码

有点奇怪的是 fwin3 的曲线在窗边缘似乎有点什么问题。