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

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

  [复制链接]
发表于 2023-8-28 12:46:33 | 显示全部楼层
本帖最后由 Jack315 于 2023-8-28 13:47 编辑

实现搜索平均意义上拟合效果最佳区域的方法:
1. 任意给定一个搜索起点。
2. 在这个点为中心的一个区域内、大致均匀的各个点上。
分别计算代表拟合优劣的一组数值。
3. 根据这组数值建立一个被平均了的、光滑的曲面模型。
4. 用这个模型计算出梯度,这个梯度就指向
极值(最大或最小)点变化最快的方向。
5. 如果搜索区域内已经包含了最小值,迭代结束;
否则根据梯度方向,朝向最小值方向移动,继续搜索的过程。

暂且称这个搜索为全局搜索。
这个全局搜索提高了找到全局最优解的可能。
可以尝试几个不同的搜索起始点,然后再作比较和判断。

在全局搜索完成后,再用软件命令找出局部最优解。

点评

软件不是我善长的,希望先生有所突破。  发表于 2023-8-28 12:49
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-8-30 00:21:05 | 显示全部楼层
本帖最后由 笨笨 于 2023-8-30 08:55 编辑
Ickiverar 发表于 2023-8-18 06:23
对于你 https://bbs.emath.ac.cn/thread-19045-1-1.html 帖子中加入一个新参数之后的情况,我用 kabcd 这五 ...

先生还在吗,楼主没看明白这一步,为啥引进红线部分函数:Sin[Pi/2 x],它有什么特殊之处,而不是其它函数?原理是啥?

  1. xlist = N@Table[Sin[Pi/2 x], {x, 2/5, 1, 1/100}];
复制代码

另外请问下面这些代码,前辈是现场直接手动一个一个打出输入到Mathematica中的?
还是说在其它软件里提前编辑好,然后把转化成的代码拷贝到Mathematica中后作小浮动修改?
比如我们在Mathtype里编辑好公式,然后直接拷贝到论坛,那么公式自然就转化成latex代码了,请问编辑下面代码有啥技巧?或者有其它软件作编程辅助,求告之。

  1. grad[f_, x0_, l0_, d_] :=
  2.   Module[{n = Length[x0], x = x0, dlist, l = l0, y, oldy, oldd,
  3.     step = 0, nkernel = Max[1, Length[Kernels[]]], ymin = Infinity,
  4.     xmin},
  5.    oldy = f @@ x;
  6.    oldd = None;
  7.    While[True,
  8.     ++step;
  9.     dlist =
  10.      Table[f @@ (x + #) - f @@ (x - #) &@
  11.         Table[d Boole[i == j], {j, n}], {i, n}]/(2 d);
  12.     dlist /= Norm[dlist];
  13.    
  14.     (*x-=l dlist;
  15.     y=f@@x;*)
  16.     y = Table[{i, #, f @@ #} &[x - l*1.5^(i - 1) dlist], {i, nkernel}];
  17.     y = SortBy[y, Last][[1]];
  18.     l *= 1.5^(y[[1]] - 1 - Round[(nkernel - 1)/2]);
  19.     x = y[[2]];
  20.     y = y[[3]];
  21.    
  22.     If[y < ymin, ymin = y; xmin = x;];
  23.     If[Mod[step, 100] == 1,
  24.      Pause[0.1]; Print[{{ymin, xmin}, {y, x}, l, step}]];
  25.     If[! (oldd === None), l *= 0.85 + 0.35 dlist.oldd];
  26.     oldd = dlist;
  27.     oldy = y;
  28.     ]
  29.    ];
复制代码

点评

刚看到论坛的坛规,原来还有这一条,受教了。  发表于 2023-8-30 09:05
没必要说话这么歹毒吧,源代码楼上文件已有,楼主发表的图片不是运行该代码,是代码中运行的理论,求解说一下,何罪之有。  发表于 2023-8-30 08:53
nyy
https://bbs.emath.ac.cn/thread-5279-1-1.html  发表于 2023-8-30 08:39
nyy
我代表gxqcn恨死你了,因为你贴源代码截图  发表于 2023-8-30 08:38
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-8-30 06:51:57 | 显示全部楼层
本帖最后由 Jack315 于 2023-8-30 07:12 编辑

先说个直觉,这个补偿函数可能有个天花板。
直觉来自图形,不是证明。先用下列代码看下图……

首先是函数的定义:
  1. target[x_] := Hypergeometric2F1[-(1/2), -(1/2), 1, x^2]
  2. ramanujan[x_] := 1 + (3 x^2)/(10 + Sqrt[4 - 3 x^2])
  3. compensation[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)
  4. err[x_, a_, b_, c_, d_] := target[x] - ramanujan[x] - compensation[x, a, b, c, d]
复制代码

接下来是个函数拟合优劣的评估函数:
  1. sseMean[aVal_, bVal_, cVal_, dVal_] :=
  2. Module[{a = aVal, b = bVal, c = cVal, d = dVal, \[Delta], data, sse},
  3.   (*计算数据点时采用的步长*)
  4.   \[Delta] = 0.01;
  5.   (*产生数据点平方*)
  6.   data = Table[(err[x, a, b, c, d])^2, {x, \[Delta], 1, \[Delta]}];
  7.   (*计算 SSE*)
  8.   sse = Total[data];
  9.   (*返回 SSE 的平均值*)
  10.   sse \[Delta]
  11.   ]
复制代码

略微说明下评估拟合的原理:
如果能100%拟合,则两曲线将会重合;
不然的话,两个曲线之间所夹的面积就代表了失拟(曲线形状失配)的程度。
这个面积应该用积分来计算。但计算积分比较耗时,
所以改成离散点上误差平方和的平均值。
与计算两曲线间的最大误差功能类似。但计算速度上可能会更快一点。

下面是用这个函数画图的代码。函数前面加了负号,
使得最小值成为最大值,就是把曲面翻了个个,便于观察。

第一幅图:
  1. Plot3D[{-sseMean[a, b, 10, 20],
  2.   -sseMean[a, b, 20, 10],
  3.   -sseMean[a, b, 20, 40],
  4.   -sseMean[a, b, 40, 20]},
  5. {a, 1, 100}, {b, 1, 100},
  6. MeshFunctions -> {#3 &},
  7. BoundaryStyle -> Thick,
  8. AxesLabel -> Automatic,
  9. PlotLegends -> {"c,d=10,20", "c,d=20,10", "c,d=20,40", "c,d=40,20"}]
复制代码


第二幅图:
  1. Plot3D[{-sseMean[10, b, c, 20],
  2.   -sseMean[20, b, c, 10],
  3.   -sseMean[20, b, c, 40],
  4.   -sseMean[40, b, c, 20]},
  5. {b, 1, 100}, {c, 1, 100},
  6. MeshFunctions -> {#3 &},
  7. BoundaryStyle -> Thick,
  8. AxesLabel -> Automatic,
  9. PlotLegends -> {"a,d=10,20", "a,d=20,10", "a,d=20,40", "a,d=40,20"}]
复制代码


第三幅图:
  1. Plot3D[{-sseMean[10, 20, c, d],
  2.   -sseMean[20, 10, c, d],
  3.   -sseMean[20, 40, c, d],
  4.   -sseMean[40, 20, c, d]},
  5. {c, 1, 100}, {d, 1, 100},
  6. MeshFunctions -> {#3 &},
  7. BoundaryStyle -> Thick,
  8. AxesLabel -> Automatic,
  9. PlotLegends -> {"a,b=10,20", "a,b=20,10", "a,b=20,40", "a,b=40,20"}]
复制代码


第四幅图:
  1. Plot3D[{-sseMean[a, 10, 20, d],
  2.   -sseMean[a, 20, 10, d],
  3.   -sseMean[a, 20, 40, d],
  4.   -sseMean[a, 40, 20, d]},
  5. {a, 1, 100}, {d, 1, 100},
  6. MeshFunctions -> {#3 &},
  7. BoundaryStyle -> Thick,
  8. AxesLabel -> Automatic,
  9. PlotLegends -> {"b,c=10,20", "b,c=20,10", "b,c=20,40", "a,b=40,20"}]
复制代码


从图上可以看出参数的不同只影响最佳参数点的位置,并不影响最佳参数点的值。
这是直觉的由来。若要证这一点,或可从误差函数的级数展开式入手。

注:全局搜索的结果就是到了这个参数最佳点,但这个点上拟合效果并不好。

再放两个演示参数变化如何影响拟合误差的动画。
选的演示点就在全局最佳参数的区域。
第一动画 err ~ a, b :
  1. Animate[Plot[err[x, a, b, 3, 10],
  2.   {x, 0, 1},
  3.   GridLines -> Automatic,
  4.   PlotRange -> All,
  5.   Frame -> True],
  6. {a, 1, 100},
  7. {b, 1, 100},
  8. AnimationRunning -> False]
复制代码

第二个动画 err ~ c, d :
  1. Animate[Plot[err[x, 9, 16, c, d],
  2.   {x, 0, 1},
  3.   GridLines -> Automatic,
  4.   PlotRange -> All,
  5.   Frame -> True],
  6. {c, 1, 100},
  7. {d, 1, 100},
  8. AnimationRunning -> False]
复制代码


在动画演示的过程中还可以看到计算发生困难的情况。
这是因为在补偿函数里采用了指数函数的原因:
\[\begin{matrix}f(x)=a^x,&a>0且a\ne1\end{matrix}\]

虽然有这个直觉,但并不能否定有可以满足拟合精度要求的解。
不过即使有,找到这个点的概率可能也不太大,就是要有点运气的意思。

鉴于成功前景不佳,所以只能另寻其它方法了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-8-30 09:50:46 | 显示全部楼层
Jack315 发表于 2023-8-30 06:51
先说个直觉,这个补偿函数可能有个天花板。
直觉来自图形,不是证明。先用下列代码看下图……

可以参靠一下楼上前辈上传的这个文件包:

Er[a,b,c,d].zip (73.51 KB, 下载次数: 0)

点评

拜读过了,貌似就是做到了(接近)天花板的样子。我这里做的是平均水平的,或者说是宏观角度的,前辈做的是局部的、精细的,离要求就差一点点  发表于 2023-8-30 10:15
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-8-30 10:21:15 | 显示全部楼层
Jack315 发表于 2023-8-30 06:51
先说个直觉,这个补偿函数可能有个天花板。
直觉来自图形,不是证明。先用下列代码看下图……


可以看看52楼吗?代码来源54楼压缩包

点评

用复制或手工输入都一样的;如果是功能,就是差不多到极致了。如果 10E-8 以下就行,还有最后一个希望:在区域打概率,因为路面不平有望中奖的。  发表于 2023-8-30 11:02
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-8-30 10:51:17 | 显示全部楼层
本帖最后由 Ickiverar 于 2023-8-30 10:55 编辑
笨笨 发表于 2023-8-30 00:21
先生还在吗,楼主没看明白这一步,为啥引进红线部分函数:Sin[Pi/2 x],它有什么特殊之处,而不是其它函数 ...


因为显然Er在x接近1的时候变化非常剧烈,需要采样更密集才能捕捉上升、下降和极值点。sinx起到让等距的x在1附近采样更密的条件。

函数是直接在Mathematica里敲的。这个是以前用过的函数所以直接粘上来了。

Mathematica是一种编程语言,不是纯粹的计算器:

它里面那些内置的优化、求解函数,FindMaximum,Solve,FindRoot什么的,都有内在的缺陷。简单函数还可以凑合用用,但不保证对复杂问题可以求解,或者解的正确性。因为它是编程语言,你需要在对问题有清晰认识的情况下,自己思考一个逻辑出来,用更基本的控制语句实现一个算法,才能充分利用Mathematica求解所有问题。

点评

还是前辈的编程妙啊,  发表于 2023-8-30 11:51
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-8-30 10:53:05 | 显示全部楼层
这是前辈区域的图:
最佳区域_01.png

这是我找到的一个区域的图:
(不能上图了

上画图代码……(参考53#)
前辈的:
  1. ContourPlot[sseMean[-2.63010099114955, 11.468733192809488, c, d],
  2. {c, 0.5, 2}, {d, 0.5, 2},
  3. Contours -> Function[{min, max}, Range[min, max, 10^-11]],
  4. PlotLegends -> Automatic
  5. ]
复制代码
  1. Plot3D[sseMean[-2.63010099114955, 11.468733192809488, c, d],
  2. {c, 0.5, 2}, {d, 0.5, 2},
  3. MeshFunctions -> {#3 &},
  4. BoundaryStyle -> Thick,
  5. AxesLabel -> Automatic,
  6. PlotLegends -> {"a,b,c,d=-2.63,11.47,0.94,1.36"}]
复制代码


这是俺的:
  1. ContourPlot[sseMean[9, 16, c, d],
  2. {c, 0.1, 10}, {d, 1, 20},
  3. Contours -> Function[{min, max}, Range[min, max, 10^-11]],
  4. PlotLegends -> Automatic
  5. ]
复制代码
  1. Plot3D[sseMean[9, 16, c, d],
  2. {c, 0.5, 7.5}, {d, 1, 10},
  3. MeshFunctions -> {#3 &},
  4. BoundaryStyle -> Thick,
  5. AxesLabel -> Automatic,
  6. PlotLegends -> {"a,b,c,d=9,16,3.6,10"}]
复制代码


注意两者数量级是一样的,只是周围地势不一样。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-8-30 11:17:10 | 显示全部楼层
Jack315 发表于 2023-8-30 06:51
先说个直觉,这个补偿函数可能有个天花板。
直觉来自图形,不是证明。先用下列代码看下图……

53#
注:全局搜索的结果就是到了这个参数最佳点,但这个点上拟合效果并不好。

这就是说,需要在这个区域用前辈的方法仔细找最佳点。

55#(点评)
还有最后一个希望:在区域打概率,因为路面不平有望中奖的。

打概率的意思就是在【最佳点区域】内,随机选取 (a,b,c,d),
然后看看拟合效果最佳的里面,有没有符合要求的。

俺找到的区域地势比较险峻,中奖的可能性或许会更大一点。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-8-30 11:44:35 | 显示全部楼层
本帖最后由 笨笨 于 2023-8-30 12:41 编辑
Ickiverar 发表于 2023-8-30 10:51
因为显然Er在x接近1的时候变化非常剧烈,需要采样更密集才能捕捉上升、下降和极值点。sinx起到让等距的x ...


发贴失误……

点评

不好意思,发贴失误  发表于 2023-8-30 12:51
说得有点过……无论如何,你的回复和回复的内容之间没有任何逻辑关系,也令我比较困惑。  发表于 2023-8-30 12:38
做出这种结论属于是对计算机、浮点数、数值计算和Mathematica的Plot的基本原理都不懂……  发表于 2023-8-30 12:21
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-8-30 12:28:57 | 显示全部楼层
本帖最后由 Ickiverar 于 2023-8-30 12:30 编辑
笨笨 发表于 2023-8-30 11:44
下面的作差效果要好点:


首先Hypergeometric的定义就是那个sum式,你这个行为和 x - x 没有区别,结果当然是恒0。

而EllipticE定义的F是定义不同的代数等价形式,它们的数值是完全精确相等的,只是因为计算式不同,因为浮点精度限制有所区别,所以画出来是毛毛糙糙的10^-16级别的误差。你可以用WorkingPrecision指定更高精度再试。

你进行任何数值计算,都会受到这个至少10^-16精度的影响。我画出10^-16级别的误差图就已经足够说明它们实际上完全一致了。

Plot[Evaluate[
  FullSimplify@
    Sum[x^(2 n) (Product[1/2 - j, {j, 0, n - 1}]/n!)^2, {n, 0,
      Infinity}] - F[x]], {x, 0, 1}, WorkingPrecision -> 30]

点评

谢谢前辈批评指导,原来如此,学习了。  发表于 2023-8-30 12:47
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-11-7 17:44 , Processed in 0.038728 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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