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

Jack315

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

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

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

笨笨

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.    ];

复制代码

Jack315

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

首先是函数的定义：

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

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

鉴于成功前景不佳，所以只能另寻其它方法了。

笨笨

Jack315 发表于 2023-8-30 06:51
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……

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

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

2023-8-30 09:50 上传

点击文件名下载附件

笨笨

Jack315 发表于 2023-8-30 06:51
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……

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

Ickiverar

笨笨 发表于 2023-8-30 00:21
先生还在吗，楼主没看明白这一步，为啥引进红线部分函数：Sin[Pi/2 x]，它有什么特殊之处，而不是其它函数 ...

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

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

Mathematica是一种编程语言，不是纯粹的计算器：

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

Jack315

这是前辈区域的图：


这是我找到的一个区域的图：
（不能上图了 ）

上画图代码……(参考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"}]

复制代码

注意两者数量级是一样的，只是周围地势不一样。

Jack315

Jack315 发表于 2023-8-30 06:51
先说个直觉，这个补偿函数可能有个天花板。
直觉来自图形，不是证明。先用下列代码看下图……

53#

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

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

55#（点评）

还有最后一个希望：在区域打概率，因为路面不平有望中奖的。

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

俺找到的区域地势比较险峻，中奖的可能性或许会更大一点。

笨笨

Ickiverar 发表于 2023-8-30 10:51
因为显然Er在x接近1的时候变化非常剧烈，需要采样更密集才能捕捉上升、下降和极值点。sinx起到让等距的x ...

发贴失误……

Ickiverar

笨笨 发表于 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]