使用9个数据点进行数值微分

Dou*_*ble 4 matlab numerical-analysis numerical-methods

当我在Matlab中尝试进行数值微分时,我遇到了问题.但我的问题可能更多的是关于数值分析而不是Matlab.

我有一个包含9个数据点的数组,表示9个不同x的f(x).我需要在数字上找到f''(x).我对x和f(x)的值是

x = [2271.38,2555.30,2697.26,2768.24,2839.22,2910.20,2981.18,3123.14,3407.06]

f(x)= [577.4063,311.3341,193.0833,141.3048,95.1501,58.8130 32.4931,6.9525,0.1481]

我可以进行插值以获得平滑的曲线.我使用样条插值但是当你想要区分时,还有其他一些插值吗?

我尝试过不同的方法:

只是简单的向前,向后和中央差异商

基于小波的方法:http://www.mathworks.com/matlabcentral/fileexchange/13948-numerical-differentiation-based-on-wavelet-transforms

和衍生套件:http://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

其中没有这个令人满意.二阶导数在步长方面非常不稳定,并且导出套件中的自适应方法工作非常糟糕.也许我只是以错误的方式使用它!

任何帮助表示赞赏!

提前致谢

小智 9

我想你前几天在MATLAB Central上问了类似的问题.你没有在那里发布你的数据,所以我当时无法给出一个好的答案.

估计二阶导数是一件困难的事情.这是一个不适定的问题.微分本身就是一个噪声放大器,所以估计二阶导数是"两次"坏.这根本不是一件容易的事,当然不是很好.

使用这组点,我选择使用我的SLM工具箱估计样条模型.

x = [2271.38, 2555.30, 2697.26, 2768.24, 2839.22, 2910.20, 2981.18, 3123.14, 3407.06];
f = [577.4063, 311.3341, 193.0833, 141.3048, 95.1501, 58.8130 32.4931, 6.9511, 0.1481];
Run Code Online (Sandbox Code Playgroud)

首先,绘制数据.我可以从那个情节中学到什么?我可能会选择做出任何推论吗?

纯数据图

简单的情节告诉我,连同你的评论,我希望这个函数是一个单调递减函数.它似乎在每一端渐近线性,如双曲线段,在整个域上具有正曲率.

所以现在我将使用此信息使用我的SLM工具箱为您的数据构建模型.

slm = slmengine(x,f,'plot','on','decreasing','on','knots',20, ...
    'concaveup','on','endconditions','natural');
Run Code Online (Sandbox Code Playgroud)

slmengine旨在以曲线形状的处方形式从您那里获取信息.您会发现,通过提供此类信息,它可以极大地规范结果的形状,以符合您对过程的了解.在这里,我只是根据您的评论对曲线形状做了一些猜测.

在上面的电话中,我指示SLM:

  • 完成后绘制结果图
  • 创建x的单调递减函数
  • 使用20个等间距的结
  • 迫使曲线到处都是正二阶导数
  • 将末尾的二阶导数设为零

生成的图是一个gui本身,允许您绘制函数和数据,也可以绘制结果的导数.垂直绿线是结点.

曲线和数据

在这里,我们看到曲线拟合是您正在寻找的合理近似值.

二阶导数图怎么样?当然,SLM是一种分段立方体工具.因此,二阶导数仅是分段线性的.这是一个问题吗?您是否会要求我提供更高阶样条的工具?太糟糕了,但不,我不会.那些高阶导数的估计太差,无法要求高度平滑的结果.事实上,我对这个预测很满意.请注意,二阶导数中的毛刺是一致的.如果我使用更多的结或更少,他们仍然在那里.这是一种很好的方式来了解您看到的形状是曲线的特征,还是仅仅是结点放置的工件.

看到我对曲线形状的限制导致了一个非常合理的拟合,尽管事实上我使用了比我的数据点更多的结.SLM在估算方面没有任何问题.

二阶导数

如果我想尝试更平滑地估计二阶导数,只需使用更多的结.SLM相对较快.因此,在50节的情况下,我们得到了二阶导数曲线非常相似的结果.

二阶导数,50节

你可以在MATLAB Central上找到SLM(这里).它确实需要优化工具箱.