Пользовательские функции в gnuplot

Перевод. Оригинал здесь.
Если вы можете написать уравнение, gnuplot может рассчитать его. Имеются в виду функции, заданные в элементарной форме, например f(x)=a*x+b и не содержащие операций интегрирования/дифференцирования, требующих численных расчетов. В gnuplot встроено много базовых функций, таких как синус, косинус, гамма-функция и т. д. Вы можете построить график любой из них, а также их комбинаций.

Для пользователей gnuplot может быть не так интересна сама возможность построения графиков функций, как использование ее для аппроксимации экспериментальных данных методом наименьших квадратов. Gnuplot может проводить также нелинейную аппроксимацию. В этом разделе мы покажем, как использовать определенные пользователем функции и как применять их для аппроксимации экспериментальных данных.

Функция, о которой пойдет речь ниже, это Лоренциан с дополнительным членом 1/sqrt(x). Это уравнение имеет четыре параметра — a, b, c и d, которые задаются экспериментальными данными.

Определение этой функции похоже на аналогичные в языках Fortran или C. В примере ниже это уранение задается в следующей форме: f(x)=c/((x-a)*(x-a)+b)+d/sqrt(x). Квадрат (x-a) также может быть записан в манере Fortran: (x-a)**2. Параметры a,b,c и d произвольны.

gnuplot> a=0.25
gnuplot> b=0.02
gnuplot> c=0.05
gnuplot> d=0.1
gnuplot> f(x)=c/((x-a)*(x-a)+b)+d/sqrt(x)
gnuplot> set xrange [0:1]
gnuplot> set yrange [0:4]
gnuplot> plot f(x)

Значения функции

Так как gnupot для построения графика рассчитывает численные значения функции, эти значения можно просмотреть с помощью команды print. Значения рассчитываются с двойной точностью.

gnuplot> print f(0.25)
2.7
gnuplot> print f(0.4)
1.33458447124371
gnuplot> a=0.4
gnuplot> print f(0.4)
2.65811388300842

Для вывода результатов в табличной форме, чтобы их можно было затем обрабатывать в других программах, используется специальный терминал table. Чтобы записать результат в файл, задайте его имя с помощью команды set output:

gnuplot> set term table
Terminal type set to 'table'
gnuplot> plot f(x)
#Curve 0, 100 points
#x y type
0 0 u
0.010101 1.63972 i
0.020202 1.39031 i
0.030303 1.30688 i
   ....

0.979798 0.191506 i
0.989899 0.188622 i
1 0.185837 i

gnuplot> set output "calc.plt"
gnuplot> replot

Поиск параметров методом наименьших квадратов

Теперь аппроксисимируем нашей функцией экспериментальные данные и найдем оптимальные значения параметров a,b,c и d. Экспериментальные данные сохранены в файле «exp.dat».

 2.5000E-03 3.0420E+00 6.47E-01
 3.5000E-03 2.5700E+00 4.37E-01
 4.5000E-03 2.3020E+00 2.53E-01
   ...

 7.0000E-01 2.7420E-01 2.14E-03
 7.5000E-01 2.5680E-01 1.81E-03
 8.0000E-01 2.4630E-01 1.59E-03

Данные состоят из трех столбцов (x,y,z), где z — абсолютная погрешность определения y, то есть имеет такую же размерность. Величина, обратная к z, представляет собой вес данной точки при аппроксимации данных. Если погрешность не задана, веса всех точек данных одинаковы.

Во-первых, построим график функции и экспериментальных данных.

gnuplot> set xlabel "Energy [MeV]"
gnuplot> set ylabel "Cross Section [b]"
gnuplot> set xtics 0.1
gnuplot> set ytics 0.5
gnuplot> plot f(x) title "Lorentzian",\
> "exp.dat" using 1:2:3 title "experiment" with yerrors

Как уже говорилось, параметры a, b, c и d произвольны, однако в данном случае их значение можно приблизительно рассчитать. Параметр a — это значение пика Лоренциана по оси абсцисс, то есть примерно 0,25 на нашем графике. Параметр b соответствует квадратному корню из ширины этого пика, в нашем случае приблизительно 0,02.

В gnuplot достаточно просто сделать аппроксимацию методом наименьших квадратов. Просто используйте команду fit и добавьте параметры с помощью опции via. Если ваша функция сильно нелинейная, необходимо осторожно подходить к выбору начальных значений параметров. В данном примере мы использовали значения, найденные выше.

gnuplot> fit f(x) "exp.dat" using 1:2:3 via a,b,c,d
 
 
Iteration 0
WSSR        : 96618.1           delta(WSSR)/WSSR   : 0
delta(WSSR) : 0                 limit for stopping : 1e-05
lambda    : 1150.73
 
initial set of free parameter values

    ...

After 17 iterations the fit converged.
final sum of squares of residuals : 3341.93
rel. change during last iteration : -5.29173e-06
 
degrees of freedom (ndf) : 47
rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : 8.43237
variance of residuals (reduced chisquare) = WSSR/ndf : 71.1049
 
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
a               = 0.26191          +/- 0.005759     (2.199%)
b               = 0.00251445       +/- 0.0008358    (33.24%)
c               = 0.00541346       +/- 0.0009206    (17.01%)
d               = 0.182469         +/- 0.007329     (4.016%)
 
 
correlation matrix of the fit parameters:
 
               a      b      c      d
a               1.000
b               0.042  1.000
c              -0.229  0.783  1.000
d               0.210 -0.538 -0.768  1.000
gnuplot> replot

Как видно по графику, наша аппроксимирующая функция недостаточно точно соответствует экспериментальным данным. Для повышения точности модифицируем член d/sqrt(x) путем введения нового параметра e, тогда он примет вид d*x*e. Примем начальное значение e=-0,5.

gnuplot> e=-0.5
gnuplot> f(x)=c/((x-a)*(x-a)+b)+d*x**e
gnuplot> fit f(x) "exp.dat" using 1:2:3 via a,b,c,d,e

     ...
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
a               = 0.25029          +/- 0.002106     (0.8412%)
b               = 0.00197707       +/- 0.0002747    (13.89%)
c               = 0.00550098       +/- 0.0003662    (6.657%)
d               = 0.21537          +/- 0.003743     (1.738%)
e               = -0.358371        +/- 0.0115       (3.208%)
 
 
correlation matrix of the fit parameters:
 
               a      b      c      d      e
a               1.000
b               0.021  1.000
c              -0.078  0.788  1.000
d              -0.110 -0.384 -0.500  1.000
e              -0.304  0.198  0.335  0.381  1.000
gnuplot> replot

Все еще не идеально. Примем a=0,24 и произведем поиск только для via b, c, d,e.

Результат выведем в файл Postscript.

gnuplot> set linestyle 1 lt 1 pt 7
gnuplot> set linestyle 2 lt 2 lw 3
gnuplot> set size 0.6,0.6
gnuplot> set term postscript eps enhanced color
Terminal type set to 'postscript'
Options are 'eps enhanced color dashed defaultplex "Helvetica" 14'
gnuplot> set output "exp.ps"
gnuplot> plot"exp.dat" using 1:2:3 title "experiment" with yerrors ls 1,\
>        f(x) title "Lorentzian" with line ls 2