Подгонка Gnuplot с модифицированными функциями Бесселя

Лучшая подходящая стратегия - обычно сводить вашу нелинейную или неполиномиальную задачу к линейной или полиномиальной. В частности, линейная задача всегда имеет одно-единственное решение. Таким образом, мы идеально подошли бы f(x) = A*x + B, где B = b * besy1(b) - это для функций Бесселя второго типа, см. Ниже правку для модифицированных функций Бесселя второго типа, которые недоступны в Gnuplot. Вы делаете это так:

fit A*x + B "datafile" via A, B

Получив B, вы можете найти b, который соответствует пересечению y = x * besy1(x) с B на x = b. Поскольку besy1(x) является колебательным, у вас может быть несколько результатов, но в зависимости от диапазона, в котором указаны ваши данные, вы можете выбрать правильный. Допустим, вы получили B = 1.2 из подгонки, тогда пересечения в интервале [0:10] будут следующими:

plot [0:10] x*besy1(x), 1.2

Если ваш интересующий регион находится около x = 4.65, где есть приблизительное местоположение одного из перекрестков, ищите точное перекресток. Расстояние между x * besy1(x) и B в этой области будет приближаться к нулю, поэтому квадрат расстояния можно аппроксимировать параболой с четко определенным минимумом:

plot [4.6:4.7] (x*besy1(x)-1.2)**2

Ваш оптимальный x = b - это положение этого минимума. Вы можете экспортировать это как данные и подогнать под параболу f(x) = a2*x**2 + b2*x + c2 с минимумом, указанным в f'(x) = 0, то есть x = -b2 / (2.*a2):

set table "data_minimum"
plot [4.6:4.7] (x*besy1(x)-1.2)**2
unset table
fit [4.6:4.7] a2*x**2 + b2*x + c2 "data_minimum" via a2,b2,c2
print -b2/2./a2

Это дает x = 4.65447163370989 для положения минимума, что соответствует оптимальному b в B = b*besy1(b).

Точность этого будет зависеть от качества квадратичной аппроксимации, которая, в свою очередь, будет зависеть от того, насколько узким является ваш диапазон значений x. В этом случае диапазон [4.6:4.7] привел к тому, что квадратичное соответствие было неплохим, но не идеальным (вы можете сузить его еще больше):

plot [4.6:4.7] "data_minimum" t "data", a*x**2+b*x+c t "quadratic fit"

Изменить

Для модифицированных функций Бесселя второго типа или других сложных функций, недоступных в Gnuplot, вы можете использовать внешний синтаксический анализатор. Например, см. Мой ответ о том, как использовать внешний код Python для анализа функций: Передача Функции Python для Gnuplot.

Вы можете использовать scipy для доступа к своей функции, изменяя скрипт Python из моего другого ответа (имя файла test.py):

import sys
from scipy.special import kn as kn

n=float(sys.argv[1])
x=float(sys.argv[2])

print kn(n,x)

и в Gnuplot используйте это как

kn(n,x) = real(system(sprintf("python test.py %g %g", n, x)))

тогда вся вышеуказанная процедура работает, просто заменяя besy1(x) на kn(1,x).