2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте его в существующую тему, а создайте новую в корневом разделе "Помогите решить/разобраться (М)".

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

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

Обязательно просмотрите тему Правила данного раздела, иначе Ваша тема может быть удалена или перемещена в Карантин, а Вы так и не узнаете, почему.



Начать новую тему Ответить на тему На страницу 1, 2, 3  След.
 
 Подгонка функции (curve fitting)
Сообщение07.08.2015, 09:25 


27/03/12
23
Есть такой вопрос.
У меня имеется осциллограмма(набор точек в файле) - сигнал с резистора в цепи RLC при разрядке конденсатора.
Сигнал представляет затухающую синусоиду.
$y(t)=I_{0}\cdot e^{-\alpha \cdot t}\sin (\omega\cdot t)$
Неизвестные коэффициенты:
$I_{0}, \alpha, \omega$
Точек много, но для решения задачи можно брать не все. а, скажем, каждую сотую точку.
Пробовал решать с помощью МНК, но застрял после того, как продифференцировал функцию невязок.
(Как ее дальше решать, я не знаю. А тем более программно).

$\sigma(I_{0},\alpha,\omega)=\sum_{i=1}^{n}(y_{i}-I_{0}\cdot e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))^{2}$
Получается вот такая система уравнений.
$\frac{\partial}{\partial I_{0}}\sigma(I_{0},\alpha,\omega)=-2\sum_{i=1}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))\cdot e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i})=0$
$\frac{\partial}{\partial \alpha}\sigma(I_{0},\alpha,\omega)=2\sum_{i=1}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))\cdot I_{0}\cdot t_{i}\cdot e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i})=0$
$\frac{\partial}{\partial \omega}\sigma(I_{0},\alpha,\omega)=-2\sum_{i=1}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))\cdot I_{0}\cdot t_{i}\cdot e^{-\alpha \cdot t_{i}}\cos (\omega\cdot t_{i})=0$

Я хочу написать программу на языке программирования, а не просто решить в каком-то пакете.
В одном из математических пакетов есть NonLinear Curve Fit, там делается большое число итераций, после чего считаются коэффициенты. Вроде бы используется метод Левенберга-Марквардта.
Буду благодарен за подсказку в том, как решить поставленную выше задачу.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 10:09 
Аватара пользователя


21/01/09
3773
Дивногорск
Выпрямить и по максимальным и минимальным точкам все найти.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 10:15 


27/03/12
23
что значит "выпрямить"?
Нужно именно для данного набора точек.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 11:33 
Заслуженный участник
Аватара пользователя


11/03/08
7940
Москва
Если точек достаточно много (ну, скажем, с десяток на период), то максимум (минимум), определяемый, как $x_{i-1}\le x_i\ge x_{i+1}$ и аналогично для минимума, будет достаточно близок к максимальным/минимальным значениям функции в окрестности точки i. Тогда, составив последовательность из максимумов и (абсолютных величин) минимумов и соответствующих им моментов времени, и прологарифмировав значения минимумов и максимумов, получим линейную регрессию, в которой коэффициент при времени будет соответствовать показателю экспоненты. А разделив период времени на число минимумов и максимумов в нём, получим период колебаний (половину периода, уточняет сержант-вычислитель Равно из батареи капитана Очевилность). Другая оценка периода - через число пересечений нуля.
Это довольно грубая оценка, но может оказаться достаточно точной для практики. Грубость её связана прежде всего с тем, что отобранные так точки выборки не попадают в точности на максимумы и минимумы, хотя, если их много, лежат достаточно близко и численно близки к максимумам и минимумам функции. Кроме того, принимается, что в максимуме синус, как сомножитель функции, равен единице (минус единице в минимуме), но ввиду экспоненциального множителя это не так, хотя при небольшом показателе экспоненты это работает, как приближение.
Даже если точность такой оценки недостаточна, она может быть использована, как начальное приближение (в том же Левенберге-Марквардте)

-- 07 авг 2015, 12:17 --

Теперь по Левенбергу-Марквардту.
Начинаем с какого-то приближения, получаем оценку $y(t)$, подставив в выражение приближенные значения параметров, находим отклонения фактических значений от полученных.
Затем пользуемся линейным приближением для этого отклонения
$\Delta y=Z\Delta a+\varepsilon$
где $\Delta a$ - отклонения истинных значений параметров от приближённых, $\varepsilon$ - включает как собственно ошибки измерения и пр., так и влияние неточности линейной модели, Z - матрица из значений производных в точках.
Далее обычная линейная регрессия даёт $\Delta a=(Z^TZ)^{-1}Z^T\Delta y$, позволяя уточнить значение вектора параметров, но значения производных обычно коррелированы, и обращать придётся почти вырожденную матрицу. Поэтому к её диагонали прибавляют константу $\lambda$ (метод Левенберга) или умножают диагональный элемент на $(1+\lambda)$ (собственно метод Левенберга-Марквардта). Чем больше лямбда - тем устойчивее численно, но медленнее (при очень большой получается метод наискорейшего спуска)

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 12:35 


27/03/12
23
Как-то не очень понятно про метод Марквардта. Сложно для реализации.
Может быть, можно из той системы уравнений, что я привел выше, выразить неизвестные и решить
систему нелинейных уравнений методом простой итерации?

$I_{0}^{(k+1)}=\varphi_{1}(I_{0}^{(k)}, \alpha^{(k)}, \omega^{(k)})$
$\alpha^{(k+1)}=\varphi_{2}(I_{0}^{(k)}, \alpha^{(k)}, \omega^{(k)})$
$\omega^{(k+1)}=\varphi_{3}(I_{0}^{(k)}, \alpha^{(k)}, \omega^{(k)})$

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 14:27 


29/09/06
4552
Мои правки к тому, что Mavlik в сообщении #1043204 писал(а):
Сигнал представляет затухающую синусоиду.
$y(t)=I_{0}\cdot e^{-\alpha \cdot t}\sin (\omega\, t{\color{magenta}+\varphi})$
Неизвестные коэффициенты:
$I_{0}, \alpha, \omega, \color{magenta}\varphi$
По-моему, нельзя для таких данных игнорировать сдвиг, 4-й параметр модели.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 14:50 


27/03/12
23
Хорошо. Нельзя игнорировать сдвиг. Но это не меняет сути.
(Хотя сигнал может начинаться и с нуля).
Попробовал программу набросать - вылетает.
Видимо, метод простых итераций не подходит.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 15:12 
Аватара пользователя


21/01/09
3773
Дивногорск
Mavlik в сообщении #1043215 писал(а):
что значит "выпрямить"?
Нужно именно для данного набора точек.

Выпрямить, значит взять модуль.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 17:34 
Аватара пользователя


21/01/09
3773
Дивногорск
Алексей К. в сообщении #1043277 писал(а):
Мои правки к тому, что Mavlik в сообщении #1043204 писал(а):
Сигнал представляет затухающую синусоиду.
$y(t)=I_{0}\cdot e^{-\alpha \cdot t}\sin (\omega\, t{\color{magenta}+\varphi})$
Неизвестные коэффициенты:
$I_{0}, \alpha, \omega, \color{magenta}\varphi$
По-моему, нельзя для таких данных игнорировать сдвиг, 4-й параметр модели.

Можно начать с первого нулевого значения функции. А сдвиг по времени учесть потом при нахождении $I_0$.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 18:38 
Заслуженный участник
Аватара пользователя


11/03/08
7940
Москва
Со сдвигом - я бы рекомендовал перепараметризовать, используя тождество $A\sin(\omega t+\varphi)=B_0\sin \omega t+ B_1\cos \omega t$

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение07.08.2015, 19:30 


27/03/12
23
Можно из вышеприведенной системы выразить первое слагаемое суммы и из него взять выражение для каждой переменной. Обозначим некоторый повторяющийся множитель как
$sz$

Далее вот так

$sz=\frac{\sum_{i=2}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))}{I_{0}e^{-\alpha \cdot t_{1}}\sin (\omega\cdot t_{1})-y_{1}}$

$\alpha=\frac{1}{_{t_{1}}}\cdot \ln(\frac{sz}{sin (\omega\cdot t_{1})})$

$I_{0}=\frac{sz}{sin (\omega\cdot t_{1})\cdot t_{1}\cdot e^{-\alpha \cdot t_{1}}}$

$\omega=\frac{1}{t_{1}}\cdot \arccos\left ( \frac{sz}{sin (\omega\cdot t_{1})\cdot I_{0}\cdot t_{1}\cdot e^{-\alpha \cdot t_{1}}} \right )$

А вот программа на паскале. Но вылетает. Возможно, метод так нельзя применять.

Код:
program mnk;

uses
  SysUtils,math;

var
  xt,yt:array[1..5000] of single;
  x,dx:single;
  I0,omega,alpha:single;
  I0k,omegak,alphak:array[1..5000] of single;
  i,j,n:integer;
  sum,sum1:single;
  sz,znam:single;

begin
I0:=2.4;
alpha:=0.1;
omega:=3.3;
x:=0;
dx:=0.1;
n:=100;

for i:=1 to n do
begin
x:=x+dx;
xt[i]:=x;
yt[i]:=I0*exp(-alpha*xt[i])*sin(omega*xt[i]);
end;

I0k[1]:=10;
alphak[1]:=6.12;
omegak[1]:=1;

for j:=1 to 1000 do
begin
sum1:=0;
for i:=2 to n do
  begin
   sum1:=sum1+(yt[i]-I0k[j]*exp(-alphak[j]*xt[i])*sin(omegak[j]*xt[i]));
  end;
znam:=-(yt[1]-I0k[j]*exp(-alphak[j]*xt[1]));
sz:=sum1/znam;

alphak[j+1]:= (-1/xt[1])*ln(abs(sz/sin(omegak[j]*xt[1])));
I0k[j+1]:=sz/xt[1]/exp(-alphak[j]*xt[1]);//sin(omegak[j]*xt[1]);
omegak[j+1]:=(1/xt[1])*arccos(sz/I0k[j]/xt[1]/exp(-alphak[j]*xt[1]));

end;
end.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение08.08.2015, 13:09 
Заслуженный участник
Аватара пользователя


11/03/08
7940
Москва
0. Я бы рекомендовал всё же включить в модель фазовый сдвиг, как ещё один параметр. Даже если Вы можете гарантировать, что в начальный момент аргумент синуса равен нулю - начальный момент процесса и первое снятие его значения не обязаны совпадать.
1. Я не вижу никаких оснований ожидать того, что выписанные Вами выражения будут сходиться. Метод простой итерации сходится при довольно узких условиях.
2. Скажите, а что будет, если Вы в выражение для sz подставите точные значения параметров, а ошибка при первом измерении будет 0? И что будет, если ошибка при нём не 0?
2. А что сложного в Марквардте? Производные Вы считать можете, перемножить матрицы недорого, обратить немногим сложнее.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение08.08.2015, 18:49 


27/03/12
23
Евгений Машеров
Я не знаю, что будет при подстановке точных значений.
А метод ньютона тоже не подойдет?

Где можно почитать про метод марквардта?

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение09.08.2015, 09:00 
Заслуженный участник
Аватара пользователя


11/03/08
7940
Москва
0. Если вдруг подставленные в формулу параметры окажутся равны точным значениям, а ошибка измерения при первом наблюдении окажется нулевой - у Вас в выражении для sz будет деление на ноль. А если ошибка не ноль - то деление на величину ошибки. То есть результат закинет неизвестно куда.
1. Ньютона - подойдёт. Только если Вы не будете учитывать корреляции между значениями производных, то сходиться будет медленно и печально. Или не будет. А если учтёте - получится Левенберг-Марквардт.
2. Лично я прочёл в книге Е.З.Демиденко "Линейная и нелинейная регрессии". Но он во многих местах описан.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение11.08.2015, 15:48 


27/03/12
23
Вот, нашел описание метода.
http://optimizaciya-sapr.narod.ru/bez_mnogomer/makva.html
Получается, что нужно решать каждое уравнение системы, и для каждого будет своя матрица Гессе по всем четырем переменным?

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 32 ]  На страницу 1, 2, 3  След.

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group