Найти коэффициенты дифференциального уравнения

[Simply me]

Добрый день! Помогите, пожалуйста, решить задачу.
Есть уравнение: $y'''+a_1y''+a_2y'+a_3y=b_0u''+b_1u'+b_2u$
Начальные условия: $y^{(i)}(0)=0, i=1,2,3$
Известны $u(t_j),y(t_j)$, j=1,...,N.
Нужно найти коэффициенты $a_1, a_2, a_3, b_0, b_1, b_2$
Пробовала заменить производные разностными аппроксимациями дифференциальных операторов:
$\frac{1}{h^3}(y(t_4)-3y(t_3)+3y(t_2)-y(t_1))+a_1(\frac{1}/{h^2}(y(t_3)-2y(t_2)+y(t_1)))+a_2\frac{1}{h}(y(t_2)-y(t_1))+a_3y=b_0\frac{1}{h^2}(u(t_3)-2u(t_2)+u(t_1))+b_1\frac{1}{h}(u(t_2)-u(t_1))+b_2u$
А как дальше решать?

[zvm]

Simply me в сообщении #827154 писал(а):

Добрый день! Помогите, пожалуйста, решить задачу.

Смотрите "Метод наименьших квадратов", "Регрессионный анализ".

[Simply me]

Спасибо. Если я правильно понимаю, метод наименьших квадратов и аппроксимация позволяют найти приближенное решение. А можно найти точное решение этой задачи и каким методом?

[zvm]

Я думал у вас задача из практики. А если значения $u(t_j),y(t_j)$, $j=1,...,N$ известны точно, то задача сводится к решению системы N линейных уравнений (для $j=1, 2,...,N$) с шестью неизвестными.

[Simply me]

На самом деле у меня измеренные значения y и u состоят из реальных значений и шума. Поэтому применение метода наименьших квадратов дает большую погрешность. А можо как-то по другому решить?

[zvm]

Simply me в сообщении #830073 писал(а):

На самом деле у меня измеренные значения y и u состоят из реальных значений и шума. Поэтому применение метода наименьших квадратов дает большую погрешность. А можо как-то по другому решить?

Вот и я про тоже. Скорее всего, у вас набор результатов измерений входного сигнала $u$ и выходного сигнала $y$ в дискретные моменты времени. Производные вы оцениваете на основе результатов измерений по стандартным формулам (или у вас есть приборы для измерения скоростей изменения величин?). Здесь и появляется основная погрешность. Причем погрешность в оценке третьей производной будет очень велика.
Но если у вас полезный сигнал тонет в шуме до почти полной неразличимости, не стоит надеяться получить достоверные результаты.
Здесь может помочь сбор очень большой статистики; как известно, доверительные интервалы искомых коэффициентов уменьшаются с ростом числа наблюдений. Ну и стоит подумать об уменьшении шумов и уменьшении интервала дискретизации.
Да! А что у вас представляют собой $u$ и $y$? Вы можете активно на них воздействовать, например, сделать $u(t)$ синусоидой с заданной частотой?

[Simply me]

Спасибо! Наверное, подумаю об увеличении числа наблюдений. Так-то метод наименьших квадратов реализовала. Может, применяя его, можно получить достаточную точность, если провести больше экспериментов. А так мне эту задачу дали больше как теоретическую. Точно не знаю, но предполагаю, что приборов для измерения производных нет. u - это характеристики нагнетательных скважин, y - добывающих.

[svv]

Возможно, будет полезной книга:
Марпл-мл. Цифровой спектральный анализ и его приложения.

Там описываются методы решения подобных задач, в которых, несмотря на дискретность отсчетов, дифференциальные операторы не заменяются разностными (и это хорошо, потому что брать численно третью производную от дискретного зашумленного сигнала — это... сами понимаете).
Один из этапов решения — определение коэффициентов $\alpha_i, \beta_i$, при которых выполняется рекуррентное соотношение, в Вашем случае имеющее вид
$y_n+\alpha_1 y_{n-1}+\alpha_2 y_{n-2}+\alpha_3 y_{n-3}=\beta_0 u_n+\beta_1 u_{n-1}+\beta_2 u_{n-2}$
для всех $n$. Здесь $y_n=y(t_n)$.

(Оффтоп)

Понятно, что выполнение ДУ влечет выполнение рекуррентных формул?
Ну, например, если $y'=ky$, то $y_n=e^{k\tau}y_{n-1}$, где $\tau=t_n-t_{n-1}$ — шаг по времени.

Если сигнал $y(t)$ только приблизительно удовлетворяет ДУ, то и рекуррентные соотношения выполняются приближенно, в этом случае коэффициенты $\alpha_i, \beta_i$ подбираются так, чтобы минимизировать, например, сумму квадратов отклонений левых частей от правых.
И уже исходя из этих коэффициентов строится дифференциальное уравнение.

Считается, что этот класс методов точнее.

[Simply me]

Спасибо!

[B@R5uk]

Вообще, лучше всего стоило бы действовать так. Для функций $y(t)$ и $u(t)$ подобрать адекватную модель, например, представить их в виде линейной комбинации небольшого (много меньшего $N$, но много большего шести) числа функций, потом потребовать минимум функции

$$f(a_1, a_2, a_3, b_0, b_1, b_2)=\int\limits_{t_1}^{t_N}{\left((y'''+a_1y''+a_2y'+a_3y)-(b_0u''+b_1u'+b_2u)\right)^2dt}$$

В таком подходе оба этапа -- это самые обычные линейные МНК. Первый МНК убъёт шум в измерениях $u(t_j)$, $y(t_j)$. Второй даст искомые параметры. Но такой подход прокатит, только если модель для функций $y(t)$ и $u(t)$ будет действительно адекватной.

-- 24.02.2014, 17:36 --

Simply me, а у вас есть конкретные данные по $y(t)$ и $u(t)$ посмотреть?

[Simply me]

Нет конкретных данных.
Уравнение на самом деле n-го порядка.
$y^{(n)}+a_1y^{(n-1)}+...+a_ny=b_0u^{(n-1)}+b_1u^{(n-2)}+...+b_{n-1}u$
А когда будет много входов, то b будет матрицей, еще добавятся слагаемые с неизвестными коэффициентами.
А почему нужно 2 раза применить МНК? Чтобы сначала найти коэффициенты функции f, а потом найти $a_1, a_2,..., a_n, b_0, b_1,...,b_{n-1}$, минимизирующие интеграл квадрата разности правой и левой части?
Еще с реализацией есть проблема. Ведь должны быть неизвестны только коэффициенты. Поэтому $y^{(n)}$ нужно заменить на что-то известное.
Вот допустим $y'(t_i)=\frac{y(t_{i+1})-y(t_i)}{h}$
Для второй производной можно вывести формулу, подставив первую производную вместо функции. Получится:
$y''(t_i)=\frac{y(t_{i+1})-2y(t_i)+y(t_{i-1})}{h^2}$
Вот для третьей:
$y'''(t_i)=\frac{y(t_{i+1})-3y(t_i)+3y(t_{i-1})-y(t_{i-2})}{h^3}$
А как получить в общем случае для n-ой производной?

[_Ivana]

Simply me в сообщении #830367 писал(а):

А как получить в общем случае для n-ой производной?

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

[zvm]

Simply me в сообщении #830104 писал(а):

u - это характеристики нагнетательных скважин, y - добывающих.

Вы знаете, в таком случае, думаю, требуется искать дифференциальное уравнение с запаздывающим аргументом:
$y'''(t)+a_1y''(t)+a_2y'(t)+a_3y(t)=b_0u''(t-\tau)+b_1u'(t-\tau)+b_2u(t-\tau)$,
$\tau$ - время запаздывания.
Предложение svv мне понравилось. Действительно, зачем искать коэффициенты дифференциального уравнения, когда можно разностного (одни в другие пересчитываются).
И тогда можно использовать следующую модель:
$y_n+\alpha_1 y_{n-1}+\alpha_2 y_{n-2}+\alpha_3 y_{n-3}=\beta_0 + \beta_1 u_{n-M}+\beta_2 u_{n-M-1}+\beta_3 u_{n-M-2}$.
Величину $M$ можно оценить, построив корреляционную функцию $y$ и $u$ и посмотрев, где она достигает максимума.
А вот предложения B@R5uk я совсем не понял. Для меня нужно подробнее.

[Simply me]

Нашла коэффициенты. Теперь возник вопрос по теории управления. Как от этого уравнения перейти к уравнению системы в пространстве состояний?

[Simply me]

Т.е. в Матлаб y1=sim(mo,u1); так идентифицируем модель в пространстве состояний, где mo=idss(A*Ts+eye(n),B*Ts,C,D,zeros(length(Xi),1), Xi);
А что передать y1est=sim(?,u1)