对 海底眼 "揭开Excel平滑曲线的秘密" 一帖的改错和威力加强

灰袍法师

xzk235 发表于 2014-5-28 00:30

我直接看 海底眼 原帖给出的英文Excel，用你的测试坐标去算，得到的坐标点也跟Excel的曲线有偏差

噫，看来得另找一个 贝塞尔插值 的代码描述，然后再修改。

至于把图表拉长以后发生的偏差，这个在原帖的英文excel有做处理，代码把图表也作为参数输入自定义插值函数，函数代码会计算一次 坐标系坐标 到 Excel图表坐标 的换算，海底眼和我的附件都没做这一步换算，所以改变图标的默认比例就会有偏差。

edwin11891

感谢好贴，深受启发！

以下为使用发现的问题，与楼主商讨。
在使用中发现，平滑后的效果还是不太理想，跟Excel自己平滑的曲线有差别，细微差别不容易验证。比如以下列数据来检验（说明一下，以下数值为土力学筛分结果，X值代表颗粒粒径，Y值代表小于该粒径的百分比，因此，在X为严格单调递减的情况下，Y也是严格单调的（至少非递增），这也是后面要验证的问题之一）：

X        100        60        40        20        10        5        2        0.5        0.25
Y        100        86        86        53        31        16        10        5        2

计算后发现，有两处有问题：第一个，在X位于40-60之间（节点2-节点3之间），拟合曲线上凸严重，见附件。这是不允许的，前面说过，曲线应该是严格单调的，这个现象，利用楼主的方法、Excel自身的平滑方法、三次样条插值，均出类似的问题。可是在 Grapher （一个优秀的绘图软件）中，利用 Spline smoothing方法进行拟合时，完美解决这个问题。可见，问题不是不能解决，还是方法有问题。可是Grapher的这个Spline smoothing方法，怎么搜都没有介绍原理。

    第二个问题，在X位于0.25-0.5之间（节点8-节点9之间），利用楼主计算方法，连接线明显偏差，Excel自身的平滑和Grapher的Spline smoothing方法均没有这个问题。这还是说明了贝塞尔方法有欠缺，需要加强改进。

     第三个问题是，利用楼主的方法将上述X、Y进行插值时，比如以0.25为步长进行插值，会有很多点得不到有效 Y 值（见附件Excel文件内 N45-O518区块），这可能是因为楼主方法内解三次方程得到的 t 值均大于1或小于0，得不到有效 t 值的缘故，但是这明显不符合事实。因为只要 X 值介于 节点1至节点9的范围内，都应该会有结果。

       由此看来，Excel自身的平滑还是独门秘技，并非贝塞尔曲线那么简单。Excel的平滑方法比贝塞尔方法好，Grapher的Spline smoothing 拟合方法比Excel更好。但是该软件中也不介绍具体的原理和方程。http://pages.cs.wisc.edu/~deboor/pgs/ 对其原理有介绍（主要为其中 smooth.f, chol1d.f, setupq.f 三个程序），Fortran功底不够，改造起来太费力。

      以上内容，与楼主参考，并再次致谢！

lee1892

edwin11891 发表于 2014-9-10 09:17
感谢好贴，深受启发！

以下为使用发现的问题，与楼主商讨。

该作者的A Practical Guide to Splines：
http://vdisk.weibo.com/s/ueXVFnVRoUDex

Fortran还是很容易看明白的，和VB很接近

灰袍法师

edwin11891 发表于 2014-9-10 09:17
感谢好贴，深受启发！

以下为使用发现的问题，与楼主商讨。

得不到有效插值的原因是 VBA计算产生的浮点误差
把 模块里面的声明语句 Const MaxEPS = 0.000000000000008
删掉一个零，变成 Const MaxEPS = 0.00000000000008 即可

至于插值跟Excel自身的平滑曲线有差异，是因为坐标轴使用了对数刻度，如果不用对数刻度，那么 贝塞尔插值 线跟Excel平滑曲线是重合的。

本来 Excel 就应该提供这个插值函数，唉。。。

lee1892

灰袍法师 发表于 2014-9-10 13:49
得不到有效插值的原因是 VBA计算产生的浮点误差
把 模块里面的声明语句 Const MaxEPS = 0.0000000000000 ...

法师，我前面贴的那本书里好像没有用贝塞尔插值，是不是里头有更好的方法计算样条曲线啊~

灰袍法师

lee1892 发表于 2014-9-10 14:57
法师，我前面贴的那本书里好像没有用贝塞尔插值，是不是里头有更好的方法计算样条曲线啊~

贝塞尔曲线是 样条曲线 的其中一种。

你贴的那本书应该是没有特意介绍 贝塞尔函数，而是介绍整类样条曲线吧。

我也没有研究过其他样条函数的算法。

edwin11891

灰袍法师 发表于 2014-9-10 17:13
贝塞尔曲线是 样条曲线 的其中一种。

你贴的那本书应该是没有特意介绍 贝塞尔函数，而是介绍整类样条曲 ...

三次样条法插值也会出现问题一，甚至节点2和节点3之间曲线上凸更严重，因此也不大好，这个我试过。

至于曲线偏离是因为使用了对数刻度，这一点倒不象合理的理由。问题是 grapher 的样条平滑方法却可以避免以上三类问题，说明还是有更好的方法。

顺便感谢 lee1892 提供了该书，前面提及的三个Fortran程序我也曾经全部改写为VBA了，但是运行总出问题，没来得及深究，也没有精力深究，因此在此求教高手。

edwin11891

灰袍法师 发表于 2014-9-10 13:49
得不到有效插值的原因是 VBA计算产生的浮点误差
把 模块里面的声明语句 Const MaxEPS = 0.0000000000000 ...

实际上，我设置为 MaxEPS = 0.000001时也还是有这样的现象出现。

另外，我将解一元三次方程的方法用盛金公式来解时，对比二者，总会有少数点得到的 X1 、X2、X3值会有所不同（当然绝大多数是相同的）。

lee1892

edwin11891 发表于 2014-9-10 18:40
三次样条法插值也会出现问题一，甚至节点2和节点3之间曲线上凸更严重，因此也不大好，这个我试过。

至 ...

我不认为你提到的那个Smooth函数是你想要的，我直接跳到该书的相关部分，大致浏览了一下，貌似它的一个参数S是用户自行选择的，不同的s值会导致样条拟合差异十分大，你可以去看看该书中的几个不同参数s得出的结果图对比。
另外，你这样的要求估计需要把那本书上的内容大致吃透，你自己都没精力，指望别人有吗？

灰袍法师

edwin11891 发表于 2014-9-10 19:02
实际上，我设置为 MaxEPS = 0.000001时也还是有这样的现象出现。

另外，我将解一元三次方程的方法用盛 ...

把子过程 Newton_Solve 的语句
    Else
        Newton_Solve = NoRoot
        Exit Function
    End If
改为
    Else
        Newton_Solve = x0
        Exit Function
    End If
就可以解决找不到解的问题。