Бесполезность методов численного интегрирования повышенного порядка

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


Вот пример, функций задана графиком.



Метод первый, прямоугольники.


     n
I = sum   { x_i h }
    i = 0

Считаем, первое отличное от нуля значение есть при x = 4, начиная с него. Шаг равен единице. I = 1 + 3 + 2 + 0 = 6.

Метод второй, трапеции.


     n
I = sum   { 1/2 (x_{i-1} + x_i) h }
    i = 1

I = 1/2 + 4/2 + 5/2 + 2/2 + 0/2 = (1 + 4 + 5 + 2) / 2 = 6.

Метод третий, Симпсона.


     n
I = sum   { 1/6 (x_{i-2} + 4 x_{i-1} + x_i) h }
    i = 2

I = 1/6 + (4+3)/6 + (1+12+2)/6 + (3+8)/6 + 2/6 + 0/6 = (1 + 7 + 15 + 11 + 2) / 6 = 6.

Если посмотреть на общий случай, так же становиться ясно, что эти методы только добавляют фильтр (не разбираюсь в DSP, кажется это называется FIR/КИХ), но не изменяют точности интегрирования. Все, что вышло за пределы окна фильтра равно сумме измерений x_i. То есть структура такая.

+----------+      +------------+      
| Сумматор | <--- | КИХ фильтр | <--- источник данных
+----------+      +------------+

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

Итог, разницы в точности нет.

Есть ещё интересные вопросы о том как влияет фильтрация вносимая аналоговой частью (хотя бы не мгновенность выборки АЦП) системы на точность интегрирования. Но это тему пока оставлю т.к. сам никогда сильно не копал.

Скатился я глубоко в теорию, в следующей серии наверно будут кватернионы, если не будет возражений.
  • +7
  • 30 марта 2012, 01:29
  • amaora

Комментарии (31)

RSS свернуть / развернуть
Чем-то вы мне напомнили анекдот:
Инженер подозревает, что все нечетные числа простые. Во всяком случае, 1 можно рассматривать как простое число, доказывает он. Затем идут 3, 5 и 7, все, несомненно, простые. Затем идет 9 – досадный случай; по-видимому, 9 не является простым числом, но 11 и 13, конечно, простые. Возвратимся к 9, – говорит он, – я заключаю, что 9 должно быть ошибкой эксперимента.
+3
всё правильно заметили, но ИМХО так получается когда функция разбита на слишком малое количество дискретов. В случае, к примеру, 10битного АЦП сигнал величиной в 3 бита — это скорее всего шум. Или по-другому, измерять сигнал с дискретизацией 3 бита для цифровой обработки сигнала — это как-то не хорошо )
0
с квантованием 3 бита, может быть?
0
А откуда такая формула для метода Симпсона?
В Википедии, например, она немного другая. Там берутся двухинтервальные (трёхточечные) куски.

В чём нарисованы графики?
0
Н-да. Если посчитать интеграл по википедивски два раза, со сдвигом в одно измерение. Потом сложить и поделить на 2. (Усреднить) То результат будет равен интегралу посчитанному по формуле в этой блог-статье. Похоже это тоже было всем очевидно кроме меня. :)
0
Графики в inkscape.
0
Плюсанул.
Кватернионы тоже приветствуются.
Может переместите свои статьи в блог «Теория, измерения и расчет»? Там им самое место.
0
есть известная теорема Котельникова про соотношение частоты дискретизации и полосы сигнала для его ПОЛНОГО (то есть математически одинакового) восстановления по дискретным отсчётам.
Из дискретных осчётов сигнал восстанавливается заменой каждого отсчёта на функцию sin(x)/x с соответствующей амплитудой и потом всё суммируется.
Соответственно если теперь надо получить интеграл от сигнала, то надо все эти sin(x)/x проинтегрировать, что будет математически равно просто сумме отсчётов. Так что при подсчёте интегралов сигнала (правильного оцифрованного по Котельникову) от интерполяции никакого толка действительно нет.

Если посчитать Фурье преобразование сигнала, то значение спектра на нулевой частоте и есть интеграл от сигнала по времени. это же верно и для дискретного преобразования.
При этом какие бы фильтры низких частот не накладывались на сигнал, значение интеграла они никак поменять не могут в силу того что у них усиление(затухание) на нулевой частоте == 1.
+1
  • avatar
  • _pv
  • 30 марта 2012, 13:46
Я сегодня на перерыве между лекцией по матану и лекцией по матану, решил отвлечься и почитать we. И попал на эту статью. Отвлекся, блин :)
+2
Это, хотя и не правильное доказательство, но в то же время правильное желание подхода. Автор топика, посмотрите, что такое математическая индукция, вам понравится.
0
Это пример, а доказательство я оставил желающим.
0
Небольшое уточнение (просто недавно наступил на грабли по собственной глупости).

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

С шагом по Х — это, как правило, само собой разумеется — обычно интегрируют с фиксированным шагом. А вот с шагом по У лично я сразу не понял, что там тоже должна быть «сетка».

У меня была система диф. уравнений (относительно скорости), которая решалась численно методом Рунге — Кутты и на каждой итерации алгоритма я интегрированием считал координаты. А интегрировал методом прямоугольников (прочитав данную статью).

Потом (в процессе экспериментов) решил использовать метод трапеций. И получил повышение точности. Сначала удивился, потом понял, что моя реализация Рунге — Кутты возвращает решение для данного шага как значение с плавающей точкой, а ТС говорил о фиксированной сетке.

Наверное звучит сумбурно, но, надеюсь, «кому нужно, тот поймет» и не наступит на мои грабли.
0
Итог, разницы в точности нет.
Это очень формальный подход!
Итак, пусть есть данные с АЦП (регулярная сетка значений функции), которые надо проинтегрировать. Сигнал надо сначала отфильтровать (в идеале — с получением float — другая «регулярная сетка»), а потом проинтегрировать. Вычислительные затраты на метод Симпсона не так велики. К чему эти труднопроверяемые предположения
… а начальные и конечные значения либо не важны… либо близки к нулю.
?
И ещё интересно, по мнению автора — float64 это “регулярная сетка”?
0
Нет, скорее дело в том, что начало или конец интервала интегрирования не в нуле. Другими словами, погрешность накопленная на разгоне будет скомпенсирована при торможении.
0
Упс, Вы правы, туплю. Действительно, начальное значение одной из координат не нулевое.
0
Условие на самом деле менее строгое, наверно, достаточно будет равенства начального и конечного значений.
0
«Нет, ребята демократы — только чай!»
На самом деле очень существенная ошибка. Я не понимаю Вашего сленга — разгон — торможение, но факт в том, что приведённые ТС примеры — только для линейных функций!
0
Соответственно — там метод трапеций рулит — а метод Симпсона, соответственно поддерживает! e_mc2 — ты то как на такую детскую разводку попал?!
0
Посмотрите вот эту ветку комментариев. Там есть выкладки коллеги amaora относительно его утверждений. Я в них ошибок не нашел.

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

Признаюсь, я поначалу тоже сильно скептически отнесся к данной статье (ибо всем известно с института о том, что методы высшего порядка рулят всегда) но теперь я не уверен, что «всегда». Ни в его аналитике, ни в членном моделировании (правда я моделировал расчеты в целых числах) я опровержения не нашел.
0
Я одну все же нашел. В формуле должен быть минус, а не плюс.
Алсо из нее можно сделать вывод, что результаты идентичны не только при x1=xn=0, но и при x1=-xn.
0
Только вот как заглянуть в будущее? Чтобы это условие выполнялось? А всего то нужно сдвинуть вправо первую и последнюю выборку — и никакой головной боли! Короче, я на 100% согласен с Perfer в этом вопросе.
0
Только вот как заглянуть в будущее?
Иногда будущее очевидно :) Например, мне доводилось заниматься моделированием всяких там физических процессов. Там, зачастую, (если решать задачу в лоб) сначала получают моментальное ускорение (через сумму всех действующий сил) потом интегрированием переходят к скорости, потом к координатам. И для множества задач начальная скорость равна нулю, и конечная равна тоже нулю (ибо сила трения в том или ином виде присутствует).

Я тогда применял трапеции, как компромиссный вариант по точности (меня так учили в моем кулинарном техникуме ). Кажется, что сложность вычислений между квадратами и трапециями копеечная, но если помножить это все на количество частиц (для которых мы делаем расчеты) и ΔT – то, в моем случае, эта оптимизация экономила бы часы симуляции.
0
между квадратами и трапециями
… прямоугольниками и трапециями…
0
Ну вопрос-то элементарный! Из первого условия ТС —
регулярная сетка значений функции заранее задана
сразу вытекает, что использование метода интегрирования выше первого никакого выигрыша не даст!
А из второго (начальное и конечное значение функции ==0) вытекает совпадение метода прямоугольников и трапеций. Как верно заметил Vga, это условие можно расширить: f(0)=-f(N)
… в моем случае, эта оптимизация экономила бы часы симуляции.
Вы это серьёзно? :-) Можно пример, желательно на Ассемблере?
0
Кстати, существует более удобная формула для метода трапеций, в отличие от приведённой ТС:

Видно, что если арифметически сдвинуть вправо первое и последнее значение, то получим полное совпадение с формулой прямоугольников.
Экономия «часов симуляции» при этом — для меня загадка! :-)
0
Кстати, существует более удобная формула для метода трапеций

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

Одна симуляция длилась сорок пять минут. Затупил программу — можно сорок пять минут официально нифига не делать. Посмотрел результаты, поменял что-то в коде или начальных условиях и опять 45 минут перерыв. И так изо дня в день. Мне сейчас уже тяжело приникнуть сколько там можно было бы сэкономить избавившись от дополнительного сложения и деления чисел с плавающей точкой, но думаю время одной симуляции уменьшилось бы ощутимо, а в обшей вылилось бы все в часы.
+1
Одна симуляция длилась сорок пять минут. Затупил программу — можно сорок пять минут официально нифига не делать.
Так это ж плюс!

Когда работал в одном КБ, так же шутили про моделиста-матлабщика — мол запустил симуляцию и официально пинаешь хуи)

А вообще оптимизация подобных циклов — задача достаточно важная, да.
0
Так это ж плюс!

:)
Смотря как к этому относиться. Мне нужен был результат а не процесс, поэтому это жутко напрягало. Я назвал сою работу «Решение задачи Коши для Золушки, вахтенным методом сутки-через-трое». Зато в этом есть свои плюсы — это отучает от привычки бездумно запускать код: приходилось его скрупулезно вычитывать перед каждым запуском.
0
Одна симуляция длилась сорок пять минут.
Это разве долго? Вот процедура генерации ОС RSX 11, к примеру, занимала практически весь день, что-то порядка 6-7 часов.
0
Это разве долго?
Ну, все относительно. Для меня было долго.

Для человека, в обязанности которого входила генерация образов ОС RSX 11 — мои проблемы, подозреваю, выглядели бы мелкими :)
0
Это еще смотря как много раз это дело запускать надо.
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.