用gnuplot图计算数值导数

Question

我一直在尝试用gnu图数值计算导数，使用另一个讨论中的脚本，甚至使用相同的数据文件。然而，我一直收到这样的错误：

gnuplot> d(y) = ($0 == 0) ? (y1 = y, 1/0) : (y2 = y1, y1 = y, y1-y2)
                                ^
         "prova.g", line 7: ')' expected

我不知道在这里该怎么办。有什么帮助吗？

Answer 1

这里是我收集的数值导数的一个例子。需要gnuplot >=5.0并使用文件而不是数据库锁(可能还需要对gnuplot>=4.4.0进行一些调整)。

脚本：(与gnuplot>=5.0.0合作，2015年1月)

### numerical derivatives
reset session

# create some data
MyFunction = "sin(x)/x"
set table $Data
    set samples 150
    plot [-10:10] '+' u 1:(@MyFunction) w table
unset table

DerivX(colX) = (x0=x1,x1=column(colX),(x0+x1)/2.)
DerivY(colY) = (y0=y1,y1=column(colY),(y1-y0)/(x1-x0))

set table $Deriv1
    plot x1=y1=NaN $Data u (DerivX(1)):(DerivY(2)) w table
unset table

set table $Deriv2
    plot x1=y1=NaN $Deriv1 u (DerivX(1)):(DerivY(2)) w table
unset table

set table $Deriv3
    plot x1=y1=NaN $Deriv2 u (DerivX(1)):(DerivY(2)) w table
unset table

plot $Data   u 1:2 w l lc rgb "red"       ti MyFunction, \
     $Deriv1 u 1:2 w l lc rgb "web-green" ti "1st Derivative", \
     $Deriv2 u 1:2 w l lc rgb "blue"      ti "2nd Derivative", \
     $Deriv3 u 1:2 w l lc rgb "magenta"   ti "3rd Derivative"
### end of script

结果：

​

​

添加：(版本为gnuket4.2.6，9月1日)。(2009年)

Gnuket4.2.6没有数据库锁和串行计算，但是如果没有这些特性，这里有一个麻烦的解决方案。

1. 为了举例说明，脚本创建了一个数据文件SO68198576.dat (您已经有了输入文件)
2. 将数据文件绘制为另一个文件，TEMP1跳过第一行数据
3. 使用系统命令TEMP2将文件逐行合并到另一个文件paste中(无论是在系统上的Linux上，还是在必须安装的paste上，例如从GnuWin安装的CoreUtils )。
4. 现在您可以计算分别来自第1列和第4列以及第2和第5列的两个连续数据点之间的dx和dy。
5. 由于文件有不同的长度，所以应该跳过最后一行。这可以通过系统命令head -n -2来完成。

这就是TEMP2的样子：

#Curve 0 of 1, 150 points   #Curve 0 of 1, 150 points
#x y type   #x y type
-10 -0.0544021  i   -9.86577 -0.0432646  i
-9.86577 -0.0432646  i  -9.73154 -0.0310307  i
-9.73154 -0.0310307  i  -9.59732 -0.0178886  i
...

脚本：(使用gnuket4.2.6，需要系统命令paste和head)

### numerical derivative for gnuplot 4.2.6
reset

FILE = "SO68198576.dat"
set table FILE
    set samples 150
    plot sin(x)/x
unset table

TEMP1 = FILE.'1'
TEMP2 = FILE.'2'
set table TEMP1
    plot FILE u 1:2 every ::1
unset table

system(sprintf('paste %s %s > %s', FILE, TEMP1, TEMP2))
system(sprintf('head -n -2 %s > %s',TEMP2, TEMP1))

x0(col) = (column(col)+column(col+3))/2.
dx(col) = (column(col+3)-column(col))/2.
dy(col) = (column(col+3)-column(col))/2.

plot FILE u 1:2 w lp pt 7 title "Data", \
     TEMP1 u (x0(1)):(dy(2)/dx(1)) w lp pt 7 title "1st Derivative"
### end of script

结果：(截图4.2.6)

​

​