找回密码
 欢迎注册
楼主: 笨笨

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

  [复制链接]
 楼主| 发表于 2023-9-11 05:55:09 | 显示全部楼层
本帖最后由 笨笨 于 2023-9-11 05:58 编辑
Jack315 发表于 2023-9-10 02:55
差不多想到一个方向上了……

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


要考虑到定义域:0≤x≤1
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-11 06:04:30 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-11 06:26 编辑

几种窗函数:
\[\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 的曲线在窗边缘似乎有点什么问题。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-11 07:11:21 | 显示全部楼层
【拟合迭代过程】
函数定义:
  1. f0[x_] := Hypergeometric2F1[-(1/2), -(1/2), 1, x^2]
  2. f1[x_, {a_, b_, c_}] := a + b x + c x^2
复制代码

初次拟合:
  1. f1Param = Solve[Table[
  2.     Table[SeriesCoefficient[f0[x], {x, 0.5, n}], {n, 0, 2}] ==
  3.      Table[SeriesCoefficient[f1[x, {a, b, c}], {x, 0.5, n}], {n, 0, 2}]
  4.     , {i, 3}], {a, b, c}][[1]]
复制代码

修正拟合使定义域两端误差为零:
  1. f1k = ((f0[1] - f0[0])/(f1[1, {a, b, c}] - f1[0, {a, b, c}])) /. f1Param
  2. f1b = Mean[{f0[0] - f1k f1[0, {a, b, c}], f0[0] - f1k f1[0, {a, b, c}]}] /. f1Param
  3. f1[x_] := Table[x^(i - 1), {i, 3}].Table[f1Param[[i]][[2]], {i, 3}] f1k + f1b //Simplify
复制代码

误差函数:
  1. err1[x_] := f0[x] - f1[x]
  2. Plot[err1[x], {x, 0, 1},
  3. PlotStyle -> Red,
  4. Filling -> Axis,
  5. GridLines -> Automatic,
  6. Frame -> True,
  7. PlotLabel -> "f1(x) 的拟合误差"]
复制代码

误差函数的三个零点:
  1. sol = FindRoot[err1[x], {x, 0.5}]
  2. x0s = {0, sol[[1]][[2]], 1}
  3. Plot[err1[x], {x, 0, 1},
  4. PlotStyle -> Red,
  5. Filling -> Axis,
  6. Epilog -> {PointSize[0.03], Table[Point[{x0s[[i]], 0}], {i, Length[x0s]}]},
  7. GridLines -> Automatic,
  8. Frame -> True,
  9. PlotLabel -> "误差函数 err1(x) 的零点"]
复制代码

零点间的极值点:
  1. xms = {x /. FindRoot[D[err1[x], x], {x, Mean[x0s[[1 ;; 2]]]}],
  2.   x /. FindRoot[D[err1[x], x], {x, Mean[x0s[[2 ;; 3]]]}]}
  3. yms = err1[x] /. x -> xms
  4. Plot[err1[x], {x, 0, 1},
  5. PlotStyle -> Red,
  6. Filling -> Axis,
  7. Epilog -> {PointSize[0.03],
  8.    Table[Point[{xms[[i]], yms[[i]]}], {i, Length[xms]}]},
  9. GridLines -> Automatic,
  10. Frame -> True,
  11. PlotLabel -> "误差函数 err1(x) 的极值点"]
复制代码

拟合左侧的峰:
  1. f2[x_, {a_, b_, c_}] := a + b x + c x^2
  2. f2Param = Solve[Table[
  3.     Table[SeriesCoefficient[err1[x], {x, xms[[1]], n}], {n, 0, 2}] ==
  4.      Table[SeriesCoefficient[f2[x, {a, b, c}], {x, xms[[1]], n}], {n, 0, 2}]
  5.     , {i, 3}], {a, b, c}][[1]]
  6. Plot[{err1[x], f2[x, {a, b, c}] /. f2Param}, {x, 0, 1},
  7. Filling -> {1 -> {2}},
  8. GridLines -> Automatic,
  9. PlotLegends -> "Expressions",
  10. Frame -> True,
  11. PlotLabel -> "拟合函数 f2(x)"]
复制代码

窗函数:
  1. fwin[x_, {m_, s_}] := E^-((x - m)/s)^2
  2. (*计算窗函数参数 s 的公式,w 是相邻零点的宽度*)
  3. winS[w_] := w/3
复制代码

接下来是通过加窗函数逐个拟合误差函数的极值点……
代码改乱了,上个图看看加窗拟合的效果:
加窗拟合效果.png

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-11 07:31:15 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-11 07:54 编辑

迭代过程每次都拟合最大极值点的峰。
拟合的方法与 103# 相比作了改进:
保证极值点成为新的零点的条件下,最小化区间 [0,1] 两端点的误差平方和。
这样不至于因为不断的迭代而使这两个点的误差增加。

如果最终最大误差出现在区间端点,
则用一直线进行拟合,就可以把两个端点的误差重新拉回到零。
在迭代过程中这两点其实是伪零点,
将其作为零点处理,可以方便代码的设计。

这就是用加窗实现分段拟合的方法,
因此可以达到任意想要的拟合精度。
只是拟合函数的项数会特别多(有限项),
不似 LZ 的 \(\varphi(x)\) 函数那般精巧,有点“暴(xi)力(tong)”的感觉。

如果窗函数更方,而不再是这般油头滑脚的,
相信拟合会更精确、迭代能更快地达到所需的精度,
即更短的拟合函数项。

进行迭代的代码设计差不多了,还需要点时间进行调试。
另外也要使用更方的窗函数……希望这次能成功了。

点评

找到啥拟合公式?能否给出一个具体实例来,误差10^-8一下  发表于 2023-9-11 20:43
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-9-11 20:44:47 | 显示全部楼层
elim 发表于 2023-9-9 23:28
还没适应本论坛的 MathJax 设置。所以贴个图:

然后呢?没下文啦???

点评

好的  发表于 2023-9-11 23:05
调试代码中,耐心等待结果。  发表于 2023-9-11 22:57
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-12 00:12:03 | 显示全部楼层
笨笨 发表于 2023-9-11 20:44
然后呢?没下文啦???

笨笨可不是一般的笨,而是刻意不动脑子不学习不动手的笨.

@Jack315 :  建议用切比雪夫理论一举得到支集区间上凸/凹函数\(\tau/\|\tau\|_{\infty}\)的最隹多项式一致通近(\(\|\;\|_{\infty}\le 10^{-5}\)).
这样综合起来的误差就控制在e-12.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-12 01:30:47 | 显示全部楼层
elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨,而是刻意不动脑子不学习不动手的笨.

@Jack315 :  建议用切比雪夫理论一举得到支集 ...

感谢前辈指导。切比雪夫理论的思路挺不错的,回头试试。
笨笨笨得可爱,这个题搞明白了啥意思其实不难,正好练习下编程。
笨笨其实不笨,光看这补偿函数就不简单,只是思路不一样,会有天花板。
赠人玫瑰,手留余香

点评

楼主脑子不好使,数学素养及能力有限,为事实!希望希望先生能得到很好的拟合函数来!  发表于 2023-9-13 22:04
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-9-13 22:01:31 | 显示全部楼层
elim 发表于 2023-9-12 00:12
笨笨可不是一般的笨,而是刻意不动脑子不学习不动手的笨.

@Jack315 :  建议用切比雪夫理论一举得到支集 ...

前辈你好,按照这个思路走,可否给个具体实例来????
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-9-13 22:07:58 | 显示全部楼层
我就纳闷了,对拉马努金公式再拟合,一个函数族且误差达到10^-12就不能搞定???
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-14 09:51:27 | 显示全部楼层
本帖最后由 Jack315 于 2023-9-14 09:55 编辑
笨笨 发表于 2023-9-13 22:01
前辈你好,按照这个思路走,可否给个具体实例来????


拟合函数:\(\lambda_0+\lambda_1x+\lambda_2x^2\)
窗函数为:\(\frac{1}{2tan^{-1}(\frac{b-a}{2}k)}\big(tan^{-1}[k(x-a)]-tan^{-1}[k(x-b)]\big)\)

代码已初步实现,有两个方面的难点:
1. 精确计算误差函数的零点和极值点,难在零点的计算上。
2. 拟合中相关数据较多,需要进行数据组织。

这是最新的调试状态,设置的精度要求为 1E-4,纯手工迭代:
加窗拟合调试.png
蓝色点为零点;红色点为极值点;
青色点为待拟合点;黄色点为加窗边界点。
注:代码中,定义域端点 [0,1] 同时作为零点和极值点处理。
这是编程的需要,不是数学上的含义。

从图上可以看到漏了一个零点,加窗出错了。

要不要我把代码发上来,让高人来调试?
俺是超级业余的,数学和编程都不专业,就是爱玩这个东东。

点评

辛苦了,谢谢  发表于 2023-9-14 20:44
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-12-22 14:27 , Processed in 0.045916 second(s), 24 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表