Метод наименьших квадратов и сложная функция.

[Hitmanmix]

Уважаемые форумчане, столкнулся с проблемой. Имеется некая функция $F1$, которая зависит от трех параметров (и времени).

$\[F1(t) = k \cdot \left( {1 - \frac{{{e^{\frac{{^{ - \zeta \cdot t}}}{{T1}}}}}}{{\sqrt {1 - {\zeta ^2}} }} \cdot \sin \left[ {\sqrt {1 - {\zeta ^2}} \cdot \frac{t}{{T1}} + \arctan \left( {\frac{{\sqrt {1 - {\zeta ^2}} }}{\zeta }} \right)} \right]} \right)\]$

Экспериментально определены значения этой функции $F1 $ для ряда значений времени $t$. Задача состоит в том, чтобы отыскать значения коэффициентов

$\[T1\,\,\,\,\,k\,\,\,\,\,\zeta \]$

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

Для этого я использовал МНК и Scilab.

Вот код

function [zr]=G(c,z)
zr=z(2)-c(1)*(1 - ((%e^((-c(2)*z(1))/(c(3))))/(sqrt(1-(c(2)^2))))*sin(((sqrt(1-(c(2)^2)))*(z(1)/c(3))) + atan((sqrt(1-(c(2)^2)))/c(2)))) 
endfunction

// Исходные данные

x=[3, 6, 8, 11, 14, 18, 21, 25, 29, 32, 38, 41, 47, 52, 59, 66, 71, 77, 83, 88, 93, 100, 108, 110];

y=[5.7, 20.8, 34.7, 59.2, 86.09, 122.1, 147.3, 176.68, 199.9, 212.9, 228.6, 231.8, 231.38, 226.3, 216.4, 206.7, 201.3, 197.05, 195, 194, 195.2, 196.7, 198.6, 199];


z=[x;y];
c=[100;0.5;12]; // Начальные приближения 

[a,err]=datafit(G,z,c)

Проблема в том, что Scilab находит решение лишь в том случае, когда начальные значения (С(1), С(2) и др.) заданы очень близко к истинным. Стоит начальное значение одного из коэффициентов сделать больше или меньше истинного, Scilab либо сообщает, что проиходит деление на нуль, либо, что результат является некорректным числом.

Что Вы можете мне посоветовать? Может есть другие методы и инструменты для этого случая? Я просто не знаю, что делать и куда копать. Заранее спасибо.

[dsge]

Очевидно, что $\zeta^2$ в вашем случае стремится к единице или даже больше её по модулю (или $\zeta$ стремится к нулю), и поскольку радикал в знаменателе, Scilab ругается и не хочет делить на 0. Упростите синус, разложив на синус суммы; может быть поможет.

[Александрович]

Нужно наложить ограничения на изменение $\zeta$.

[Karan]

Тема перемещена из форума «Математика (общие вопросы)» в форум «Карантин»
по следующим причинам:

- неправильно набраны формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы);

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.

[Karan]

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»
Тема возвращена

[Евгений Машеров]

С такой убывающей экспонентой оценивать вообще трудно. В ходе решения появляются матрицы, близкие к вырожденным. На измерения ошибка наложена? Насколько она велика? Можно ли что-то из коэффициентов оценить из априорных соображений? Или, может быть, наблюдаем достаточно долго, чтобы $F_1(t)\approx k$ при $t \gg 0$? Не стоит ли трактовать $\zeta=\cos \varphi$ и $\sqrt{1-\zeta^2}=\sin \varphi$, несколько упрощая выражение? И расписать сумму под синусом?

[B@R5uk]

$\[F_1(t) = k \cdot \left( {1 - \frac{{{e^{\frac{{^{ - \zeta \cdot t}}}{{T_1}}}}}}{{\sqrt {1 - {\zeta ^2}} }} \cdot \sin \left[ {\sqrt {1 - {\zeta ^2}} \cdot \frac{t}{{T_1}} + \arctg \left( {\frac{{\sqrt {1 - {\zeta ^2}} }}{\zeta }} \right)} \right]} \right)\]$

Я так понимаю, что величина $\zeta$ тут ограничена диапазоном $(0,1)$. Тогда действительно можно положить:

$\zeta=|\cos\varphi|$

$\sqrt{1-\zeta^2}=|\sin\varphi|$

Кроме того, надо хорошо понимать, что оптимизируемая программой функция должна себя хорошо вести на всей области оптимизируемых аргументов: она должна быть непрерывной (это не критичное, но желательное требование), обязана не обращаться в бесконечность и не становиться комплексной. Кроме того, в силу особенностей численной оптимизации бесконечность не может быть экстремальной точкой: при изменении параметров оптимизации в строну бесконечности либо должна быть гарантирована бесконечная расходимость, либо должно быть бесконечно много экстремумов. Именно поэтому параметр $\zeta$ приравнивается периодичной функции по параметру $\varphi$, который теперь может меняться как угодно (разные его значения эквивалентны). Принимая это всё во внимание, включим смекалку и произведём следующие преобразования:

$\sin\left(\arctg\left(\frac{\sqrt{1-\zeta^2}}{\zeta}\right)\right)=|\sin\varphi|$

$\cos\left(\arctg\left(\frac{\sqrt{1-\zeta^2}}{\zeta}\right)\right)=|\cos\varphi|$

$$\frac{1}{\sqrt{1-\zeta^2}}\sin\left(\frac{\sqrt{1-\zeta^2}t}{T_1}+\arctg\left(\frac{\sqrt{1-\zeta^2}}{\zeta}\right)\right)=\frac{1}{|\sin\varphi|} \left(\sin\left(\frac{|\sin\varphi|t}{T_1}\right)|\cos\varphi|+\cos\left(\frac{|\sin\varphi|t}{T_1}\right)|\sin\varphi|\right)=$$

$$=\frac{|\cos\varphi|t}{T_1}\operatorname{sinc}\left(\frac{|\sin\varphi|t}{T_1}\right)+\cos\left(\frac{|\sin\varphi|t}{T_1}\right)$$

Ну и там с экспонентой перед синусом. Всё, никаких комплексных чисел и деления на ноль, лишь хорошие непрерывные функции. Теперь можете оптимизировать на здоровье. Заметьте только особенность: в последней формуле функция не синус, а синк, понимаемая в математическом смысле. Всякие среды, например, MATLAB, понимают эту функцию в физическом смысле, поэтому в программе, возможно, придётся добавить пару коэффициентов.