Помогите придумать апроксимацию

[_Ivana]

B@R5uk в сообщении #872607 писал(а):

И опять же, я получил желаемую точность

Не получили, если вспомнить про производную. А я на сплайнах получил.

B@R5uk в сообщении #872607 писал(а):

Правильно?

Правильно. И расчет этих 6 параметров делается всего за несколько операций из исходных табличных данных. Например, пара из них являются напрямую одними из табличных значений. Если испытываете затруднения с выводом формул, могу скинуть. Если вам интересно - выложите 3 таблицы на интервале $[0, 0.3]$ с шагом $0.01$ - значения функции, ее первой и второй производных, и сможете проверить те несколько строчек кода, которые рассчитывают значение сплайна. В принципе можно для каждого интервала хранить не 3 значения (функция и 2 производные в правой границе), а сразу 6 - рассчитанные коэффициенты полинома - быстрее расчет (только полином), но в 2 раза больше памяти.

[B@R5uk]

_Ivana в сообщении #872615 писал(а):

Правильно.

Спасибо, теперь понял, что вы имеете в виду.

_Ivana в сообщении #872615 писал(а):

Если испытываете затруднения с выводом формул, могу скинуть.

Спасибо. Сейчас я понимаю идею, поэтому вывести смогу сам.

_Ivana в сообщении #872615 писал(а):

а сразу 6 - рассчитанные коэффициенты полинома - быстрее расчет

Да, на мой взгляд это более рационально.

[_Ivana]

B@R5uk в сообщении #872623 писал(а):

поэтому вывести смогу сам

Выведите, сравним сложность и количество операций

B@R5uk в сообщении #872623 писал(а):

Да, на мой взгляд это более рационально.

Это зависит от того, чего вам жальче - нескольких десятков байт памяти на хранение или нескольких десятков тактов процессора на вычисления

[B@R5uk]

_Ivana в сообщении #872626 писал(а):

или нескольких десятков тактов процессора на вычисления

Помножив на миллион исполнений функции, получим 10 миллионов лишних тактов.

[_Ivana]

Тогда еще проще - $100500$ миллионов точек и интерполяция нулевого порядка - просто ближайшее значение из таблицы, и никаких вычислений

[B@R5uk]

Я тут поэкспериментировал. Думаю, сплайнами гораздо приятнее будет аппроксимировать такую функцию:

$$f\left( x \right)-1+\left( {{a}_{0}}+{{a}_{1}}x+{{a}_{2}}{{x}^{2}} \right)\exp \left( -\frac{b}{x} \right)$$

$$\begin{align} & {{a}_{0}}=0.901887608276081 \\ & {{a}_{1}}=0.264627860007153 \\ & {{a}_{2}}=-0.058890554523859 \\ & b=0.735508101699491 \end{align}$$

[B@R5uk]

На отрезке $[1,2]$ целиком удалось обойтись полиномом 7-ой степени и получить точность $7.5\cdot 10^{-8}$ для ряда и $2.5\cdot 10^{-6}$ для его производной. Коэффициенты разложения получились (от старших степеней к младшим):

   0.001043984326656
  -0.000573498066944
  -0.002211849966080
   0.011791603495936
  -0.047267312881280
   0.142018271796656
  -0.284093216907936
   0.284097632438638

Центр разложения в точке $1.5$.

Полином 6-ой степени на том же отрезке с пол пинка давал точность $3\cdot 10^{-7}$. Таким образом, увеличение степени полинома на единицу увеличило точность всего на полтора порядка, а провозиться пришлось ого-го-го, прежде чем "нащупал" оптимальное решение. Вот, думаю, может стоило разбить отрезок пополам, да аппроксимировать на каждом полиномом 5-ой степени?

[B@R5uk]

В общем, закончил я решать эту задачу. Поэкспериментировав, пришёл почти к той же идее, что предложил _Ivana, а именно, требовать от аппроксимирующего полинома совпадения значений и производных на концах отрезка аппроксимации. Это удобно, потому что полиномы на соседних отрезках оказываются автоматически сшиты до нужной степени гладкости, и при оптимизации не нужно об этом дополнительно заботиться. Отличие моего подхода заключается в том, что я считаю излишним требовать совпадения второй производной аппроксимирующего полинома и аппроксимируемой функции на концах отрезка (во всяком случае, для моей задачи это не нужно). Вместо этого лучше потребовать минимума абсолютной величины невязки на отрезке аппроксимации. То есть четыре степени свободы съедаются жёсткими связями на граничные условия, а остающиеся две-три-четыре степени свободы можно использовать, чтобы получить наилучшее приближение. В моём случае это были полиномы пятой степени (шесть коэффициентов, двое из них свободны для оптимизации).

Область аргументов аппроксимируемой функции я разбил на десять участков. Для $x\in[2,+\infty]$ я оставил одно слагаемое ряда и получил аппроксимирующую функцию $g_{10}(x)=\frac{4}{\pi}\exp(-x)$. Для $x\in[0,0.137]$ я использовал аппроксимирующую функцию $g_{0}(x)=\exp(p_0(x)-\frac{b}{x})$, где $p_0(x)$ — полином второй степени, а $b$ — константа. Отрезок $[0.137,2]$ я разбил на 8 участков и в качестве аппроксимирующей функции на каждом из них использовал полиномы пятой степени $g_{k}(x)=p_{k}(x)$, $k=\overline{1,8}$. Точки разбиения подбирал так, чтобы максимальная относительная невязка на каждом отрезке была одинакова. В результате с 10-ю участками и полиномами 5-ой степени удалось достичь точности $7\cdot 10^{-8}$. Получилась такая приятная картинка:



Абсолютная погрешность производной $10^{-5}$, но для моей задачи это не страшно, тем более, что для $x\in[2,+\infty]$ она существенно меньше.

Если кому-то будет интересно поиграться с аппроксимирующей функцией, прилагаю её скрипт для среды MATLAB:

(heating_func.m)

function [y,yd1]=heating_func(x)
x=abs(x);
y=zeros(size(x));
yd1=zeros(size(x));
xsep=[
    0.0000
    0.1370
    0.2064
    0.2874
    0.4504
    0.6000
    0.7848
    1.0300
    1.4008
    2.0000
    ];
xmap=cell(10,1);
for k=1:9
    xmap{k}=find(x>=xsep(k)&x<xsep(k+1));
end
xmap{10}=find(x>=xsep(10));

a=[
   -5.750653895630414
    3.683698706777679
   -0.982526044380308
    ];
b=0.635465797673163;
xx=x(xmap{1});
tmp=exp(polyval(a,xx)-b./xx);
y(xmap{1})=1-tmp;
yd1(xmap{1})=-tmp.*(polyval(polyder(a),xx)+b./xx.^2);
yd1(x<0.01)=0.0;

a=[
  -144.81809710607  -20.02894377021    5.21276060625    1.83719955731    0.41942048980    0.05885227953    0.00112746460   -0.00187990977
    23.41159345255   -3.68258172789   -3.77654338415   -1.01861918246   -0.21268404753   -0.01442679345    0.01327728433    0.00970515107
     2.52936486153    4.06168571317    1.65335453927    0.32966423792   -0.00488599537   -0.07103282584   -0.06203529702   -0.03874020251
    -2.09004585638   -1.20123873791   -0.17319146490    0.22457590398    0.28480740793    0.25206599491    0.18851872527    0.11624746681
    -0.34287831682   -0.59392381879   -0.74298932529   -0.71923092342   -0.62958500263   -0.51275921937   -0.37756394658   -0.23250597438
     0.98529867349    0.94921049228    0.86512112038    0.74928343895    0.63626106124    0.51372345992    0.37762427873    0.23250678472
    ];

for k=2:9
    xx=x(xmap{k})-0.5*(xsep(k)+xsep(k+1));
    y(xmap{k})=polyval(a(:,k-1),xx);
    yd1(xmap{k})=polyval(polyder(a(:,k-1)),xx);
end
y(xmap{10})=4/pi*exp(-x(xmap{10}));
yd1(xmap{10})=-4/pi*exp(-x(xmap{10}));
end

_Ivana, если вас не затруднит, то буду очень благодарен, если вы проверите эту функцию на ваших данных.

[_Ivana]

Судя по прошлому опыту, проверка излишня. Единственно что хочу посоветовать - продифференцируйте численно ваш второй график и результат сравните с третьим - для проверки.

[B@R5uk]

_Ivana в сообщении #873550 писал(а):

продифференцируйте численно ваш второй график и результат сравните с третьим - для проверки.

Ну, я это могу и аналитически сделать, как производную отношения двух функций, поскольку производная каждой функции в отношении мне известна. Получающийся график разумно выглядит.

_Ivana, большое СПАСИБО за участие и конструктивные советы.

[Yahot]

А можно файлик Excel с исходными табличными данными аппрокс. функции.
Хочу прикинуть аналитическое ее описание.