определение f''(x) для функции заданной при помощи точек

[slbel]

Добрый день!

Необходимо найти вторую и первую производную для каждой точки функции заданной при помощи точек на декартовой плоскости (с целью построения графика кривизны). Все точки функции расположены равномерно. Функция достаточно гладкая (если представить ее аналитически). Точек много - порядка E4-E6.



На графике ось X - миллисек. Этот процесс идет непрерывно по времени, и нужно обсчитывать эту кривую достаточно быстро. Точки расположены через 1 мс.

Как я решаю эту задачу сейчас:
Берем первый кусок функции в 100 точек. Нахожу ее аналитическое представление ввиде полинома третьей степени (аппроксимация методом наименьших квадратов). Далее, для точек с 25 по 74 (50 точек) нахожу кривизну во всех точках. Потом находим аналитическое представление для точек с 50 по 150. И находим кривизну в точках с 75 по 125. И так далее. Беда в том, что таким способом, мы имеем небольшие разрывы в результирующей функции кривизны через каждые 50 точек. А мне бы хотелось, чтобы в данном месте она была достаточно гладкая (хотя бы первая производная).
Хотелось бы наложить дополнительное ограничение для последующих кусков ввиде совпадения первой производной в начальной точке следующего куска с производной в последней точке предыдущего куска при определении аналитического представления.
Как это можно сделать? Или как можно находить кривизну непрерывно, не разбивая ее на кусочки?

[photon]

А в чем проблема?
У Вас есть массивы $X$ и $Y$, задающие $y(x_i)$.
Тогда производные находятся как:

$y'(x_i) =\frac{y(x_i ) - y(x_{i-1} )}{\Delta x}$

$y''(x_i) = \frac{y(x_{i+1} )- 2y(x_i )+y(x_{i-1})}{(\Delta x)^2}$

[Zai]

Если Вы успеваете посчитать вторую производную в каждой точке, то заносите ее в массив из 50 значений и проводите усреднение встроенной векторной функцией (она достаточно быстрая).

[нг]

перенесено из «Помогите Решить/разобраться»

[Руст]

photon писал(а):

А в чем проблема?
У Вас есть массивы $X$ и $Y$, задающие $y(x_i)$.
Тогда производные находятся как:

$y'(x_i) =\frac{y(x_i ) - y(x_{i-1} )}{\Delta x}$

$y''(x_i) = \frac{y(x_{i+1} )- 2y(x_i )+y(x_{i-1})}{(\Delta x)^2}$

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

[Yuri Gendelman]

photon писал(а):

...производные находятся как:
$y'(x_i) =\frac{y(x_i ) - y(x_{i-1} )}{\Delta x}$
$y''(x_i) = \frac{y(x_{i+1} )- 2y(x_i )+y(x_{i-1})}{(\Delta x)^2}$

Формулы, приведенные photon'ом - это стандартные формулы, применяемые в большинстве случаев при программировании численных алгоритмов обработки экспериментальных данных, если функция достаточно гладкая. Судя по графику, где на оси ординат отложены времена от 0 до 30000 мс, эти формулы должны быть вполне работоспособны. Эти формулы ведь в сущности те же аппроксимирующие полиномы степени 1.
Потеря точности тут может возникать при очень зашумленных данных (тогда нужны фильтры ВЧ, типа скользящего среднего).

Вот если бы на том же графике на оси ординат был диапазон от 0 до 300 мсек, то тогда, как пишет Руст, такое простое приближение не годилось бы.

А если бы на том же графике на оси ординат был диапазон от 0 до 30 мсек, то тогда, при частоте Найквиста явно больше 2000, шаг 1 мсек был бы просто слишком мал для вычислений.

[незваный гость]

Мне кажется, что мы имеем дело с измеряемым сигналом. Это ведет к нескольким следствиям: (1) ограниченность точности измерения, например, 12 битами; (2) случайные ошибки измерения. Уже первое может дать достаточно серьезные скачки производной (не говоря уж о второй производной), если вычислять по мгновенным значениям.

Именно поэтому, я думаю, slbel и применяет МНК для построения полинома 3й степени.

~~~

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

1) путь photonа: считать производную и сглаживать её. В зависимости от требований, я бы применил либо Батерворта (с временной константой от 100 мс; работает мгновенно, но имеет задержку отклика. Зато можно применять в реальном масштабе времени), либо скользящее среднее (возможно, многократное; сдвига нет, но можно делать только для обработки после, либо выдавая ответ с соответствующей задержкой).

2) Ваш путь: тут есть вилка, в зависимости от того, можно ли считать узлы равноотстоящими.

2а) Для равноотстоящих узлов всё просто (почти). Построение полинома вырождается в умножении некоторых сумм на обратную к постоянной матрице (а значит, на постоянную матрицу). Т.е., задача решается со свистом… Фокус — в пересчете сумм, служащих для построения полинома на каждом шаге (ну и накоплении ошибки). Фактически при этом вычисление становится линейным фильтром с постоянными коэффициентами.

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

~~~

В целом, всё это сводится к тому, что Вам надо копать литературу по цифровой обработке сигналов (digital signal processing, DSP).

[slbel]

Всем спасибо.

Метод предложенный Photon действительно дает скачки, ни о какой гладкости речи тут быть не может

По поводу усреднения производных по 50-100 точек. Надо попробовать, честно говоря мне такая мысль не приходила в голову.

Хотелось бы подробнее про

Цитата:

2а) Для равноотстоящих узлов всё просто (почти). Построение полинома вырождается в умножении некоторых сумм на обратную к постоянной матрице (а значит, на постоянную матрицу). Т.е., задача решается со свистом… Фокус — в пересчете сумм, служащих для построения полинома на каждом шаге (ну и накоплении ошибки). Фактически при этом вычисление становится линейным фильтром с постоянными коэффициентами.

предложенный Незваный гость.
Как этот метод называется? Если можно, нормальную ссылку в инете на описание метода.

[ГАЗ-67]

Руст писал(а):

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

Попробуйте прологарифмировать числитель и знаменатель , потом вычислите разность логарифмов и потенцируйте . Может получиться вполне приемлемое значение. Это достаточно эффективный приём при работе с очень большими и очень малыми числами . А какие проблемы возникают со сплайнами ?

[незваный гость]

ГАЗ-67 писал(а):

Попробуйте прологарифмировать числитель и знаменатель…

Для «грязных» данных это не поможет.

ГАЗ-67 писал(а):

А какие проблемы возникают со сплайнами?

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

slbel писал(а):

Как этот метод называется? Если можно, нормальную ссылку в инете на описание метода.

Метод наименьших квадратов? Бог с ним, с названием.

Давайте разберемся по существу. Мы взяли $2 n + 1$ точку, и пытаемся аппроксимировать МНК: $\min \sum\limits_{k=-n}^{n} (y_{m+k}- (a x_{m+k}^3+b x_{m+k}^2+c x_{m+k} + d))^2$. С таким же успехом мы можем сдвинуть координату на $x_m$, и получить новый полином (с другими $a, b, c, d$): $ y = a (x-x_m)^3+b (x-x_m)^2+c (x-x_m)+d$. Но! МНК станет: $\min \sum\limits_{k=-n}^{n} (y_{m+k}- (a (x_{m+k}-x_m)^3+b (x_{m+k}-x_m)^2+c (x_{m+k}-x_m) + d))^2$. А, поскольку узлы равноотстоящие, мы имеем $x_{m+k}-x_m = \delta k$, и, соответственно, $\min \sum\limits_{k=-n}^{n} (y_{m+k}- (a \delta^3 k^3+b \delta^2 k^2+c \delta k+ d))^2$. Или, обозначая $A = a \delta^3$, $B = B \delta^2$, $C = C \delta$, $D = d$, $\min \sum\limits_{k=-n}^{n} (y_{m+k}- (A k^3+B k^2+C k+ D))^2$. При этом наш интерполяционный полином станет $ y = A (\frac{x-x_m}{\delta})^3+B (\frac{x-x_m}{\delta})^2+C (\frac{x-x_m}{\delta})+D$.

Теперь дифференцируем по $A, B, C, D$ и получаем уравнения: $A \sum k^3 + B \sum k^2 + C \sum k + D \sum 1 = \sum_k y_{m+k}$, $A \sum k^4 + B \sum k^3 + C \sum k^2 + D \sum k = \sum_k k y_{m+k} $, $A \sum k^5 + B \sum k^4 + C \sum k^3 + D \sum k^2 = \sum_k k^2 y_{m+k} $, $A \sum k^6 + B \sum k^5 + C \sum k^4 + D \sum k^3 = \sum_k k^3 y_{m+k} $.

(1) Очевидно, что матрица системы постоянна, в смысле от $m$ не зависит. Более того, она легко и непринужденно точно обращается любой системой компьютерной алгебры. Не хочется вбивать формулы — тоже неплохо. Можно вычислить и обратить в программе.

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

(3) Следующий фокус позволяет резко сократить объем вычислений при переходе от $m \to m+1$. Обозначим $z_{m}^{(p)} \to \sum\limits_{k=-n}^{n}k^p y_{m+k}$.

Действительно, например, $z_{m+1}^{(0)} = \sum\limits_{k=-n}^{n}y_{m+1+k} =$ $ \sum\limits_{k=-n+1}^{n+1}y_{m+k} = $ $ \sum\limits_{k=-n}^{n}y_{m+k} +y_{m+n}-y_{m-n} = $ $ z_{m}^{(0)}+y_{m+n}-y_{m-n}$. Таким образом, мы за два сложения смогли перевычислить $z^{(0)}$.

$z_{m+1}^{(p)} = \sum\limits_{k=-n}^{n} k^p y_{m+1+k} = $ $\sum\limits_{k=-n+1}^{n+1} (k-1)^p y_{m+1+k} = $ $\sum\limits_{k=-n}^{n} (k-1)^p y_{m+1+k} + n^p y_{m+1+n} -(-n-1)^p y_{m-n} = $ $\sum\limits_{j=-0}^{p} (-1)^p C^j_p z_{m}^{(p)} + n^p y_{m+1+n} -(-n-1)^p y_{m-n} $.

“All together now”: за весьма ограниченное и независящее от $n$ количество сложений/умножений мы можем перейти от интерполяционного полинома в $m$ к полиному в $m+1$. Это дает основание считать производные при $x = x_m$ (первая $B$, вторая — $2 C$) и переходить к следующему полиному.

Бич этого подхода — накопление ошибки. Всё точно складывается и вычитается только на бумаге. Реально необходимо по крайней мере двойная точность (double).

[slbel]

Спасибо

Слишком длинный ник, у тебя Незваный гость

Есть вопрос: никак не могу сходу сообразить как считать
sum k. Вроде = -n-(n-1)-(n-2)-..+0+1+2+..+(n+1)+n = 0
(прошу прощения, не совсем понимаю как пользоваться тегом math )

Цитата:

Теперь дифференцируем по $A, B, C, D$ и получаем уравнения: $A \sum k^3 + B \sum k^2 + C \sum k + D \sum 1 = \sum_k y_{m+k}$, $A \sum k^4 + B \sum k^3 + C \sum k^2 + D \sum k = \sum_k k y_{m+k} $, $A \sum k^5 + B \sum k^4 + C \sum k^3 + D \sum k^2 = \sum_k k^2 y_{m+k} $, $A \sum k^6 + B \sum k^5 + C \sum k^4 + D \sum k^3 = \sum_k k^3 y_{m+k} $.

Попробовал вариант с аппроксимацией полиномами Чебышева с ограничениями с помощью функции
BuildChebyshevLeastSquaresConstrained
взятой на
http://alglib.sources.ru/interpolation/ ... quares.php

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

Ну и метод Незваного гостя охота попробовать. Только вот суммы по к не могу сообразить как посчитать.

[photon]

slbel, а Вы не могли бы выложить файлик с самими массивами - интересно побаловаться с аппроксимацией.

slbel писал(а):

прошу прощения, не совсем понимаю как пользоваться тегом math

А Вы почитайте тут, может MathType есть? Потренеруйтесь тут

[slbel]

Для тех кто хочет попробовать аппроксимировать, я выложил исходные данные в архиве src_approx_data.rar
на dump.ru (номер файла n8115505, см. cсылку ниже). К нему идет комментарий (см.файл comment.txt в архиве).

http://dump.ru/files/n/n8115505/

Исходные данные к графику в начале темы находятся в файле data2.txt.

[worm2]

Цитата:

2 столбец - значение по Y. Функция была предварительно отфильтрована медианным фильтром для снижения шума.
3 столбец - значения уже построенной аппроксимации по методу МНК с применением Чебышевских полиномов 5-ой степени (см. функцию BuildChebyshevLeastSquaresConstrained). Аппроксимация кусочная, удовлетворяет требованию C1.

Вам бы 3-й столбец сдвинуть вправо миллисекунд эдак на 130. Тогда совпадение со вторым будет почти идеальным.

[slbel]

Ключевое слово "почти".
Вы попробуйте построить оба графика.