Применение сплайнов для численного интегрирования

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
(Добавил основную часть теории)
(категория)
 
(7 промежуточных версий не показаны.)
Строка 1: Строка 1:
== Введение ==
== Введение ==
-
 
Ставится задача вычислить интеграл вида
Ставится задача вычислить интеграл вида
{{eqno |1}}
{{eqno |1}}
Строка 6: Строка 5:
где <tex>a</tex> и <tex>b</tex> - нижний и верхний пределы интегрирования; <tex>f(t)</tex> - непрерывная функция на отрезке <tex>[a,b]</tex>.
где <tex>a</tex> и <tex>b</tex> - нижний и верхний пределы интегрирования; <tex>f(t)</tex> - непрерывная функция на отрезке <tex>[a,b]</tex>.
-
Введем на отрезке интегрирования равномерную сетку, определим значения функции в узлах сетки. Пусть имеется совокупность узлов
+
Пусть на отрезке интегрирования дана сетка, т.е. имеется совокупность узлов
-
<tex>\left\{{t_i}\right\}_{i = 0}^{N},\; t_i = a + i{\tau},\; {\tau}= (b - a)/{N}, \;t \in \left[{a, b}\right].</tex>
+
<tex>\left\{{t_i}\right\}_{i = 0}^{N},\; i\in \left[{a, b}\right].</tex> Тогда интервал <tex>[a,b]</tex> разобьется на участки <tex>\tau_i = t_i-t_{i-1},\; i=1,\dots,N.</tex>.
-
Пусть также задана таблица
+
Пусть также известны значения функции в узлах этой сетки, т.е. задана таблица
<tex>f_i = \left\{{f(t_i)}\right\}_{i = 0}^{N}.</tex>
<tex>f_i = \left\{{f(t_i)}\right\}_{i = 0}^{N}.</tex>
Представим интеграл {{eqref|1}} в виде суммы интегралов по частичным отрезкам:
Представим интеграл {{eqref|1}} в виде суммы интегралов по частичным отрезкам:
Строка 14: Строка 13:
::<tex>\int_a^bf(t)dt=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}f(t)dt.</tex>
::<tex>\int_a^bf(t)dt=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}f(t)dt.</tex>
-
Сущность большинства методов вычисления определенных интегралов состоит в замене подынтегральной функции <tex>f(t)</tex> на отрезке <tex>[t_{i-1},\;t_i]</tex> аппроксимирующей функцией <tex>\varphi(t)</tex>, для которой можно легко записать первообразную в элементарных функциях, т.е.
+
Сущность большинства методов вычисления определённых интегралов состоит в замене подынтегральной функции <tex>f(t)</tex> на отрезке <tex>[t_{i-1},\;t_i]</tex> аппроксимирующей функцией <tex>\varphi(t)</tex>, для которой можно легко записать первообразную в элементарных функциях, т.е.
{{eqno |3}}
{{eqno |3}}
-
::<tex>\int_{t_{i-1}}^{t_i}f(t)dt=\int_{t_{i-1}}^{t_i}\varphi(t)dt+R=S+R,</tex>
+
::<tex>\int_{t_{i-1}}^{t_i}f(t)dt=\int_{t_{i-1}}^{t_i}\varphi(t)dt+R_i=S_i+R_i,</tex>
-
где S - приближеное значение интеграла; R - погрешность вычисления интеграла. Лучше всего изучена замена <tex>f(t)</tex> алгебраическим многочленом.
+
{{eqno|1*}}
 +
::<tex>J=\sum_{i=1}^N \left( \int_{t_{i-1}}^{t_i}\varphi(t)dt+R_i \right)=S+R,</tex>
 +
где <tex>S_i</tex> - приближённое значение интеграла на i-м отрезке, <tex>S=\sum_{i=1}^NS_i,</tex> а <tex>R_i</tex> - погрешность вычисления интеграла, <tex>R=\sum_{i=1}^NR_i,</tex>. Лучше всего изучена замена <tex>f(t)</tex> алгебраическим многочленом.
== Изложение метода ==
== Изложение метода ==
-
 
+
=== Общий случай ===
-
Возьмем в {{eqref|3}} в качестве аппроксимирующей функции кубический сплайн:
+
Возьмём в {{eqref|3}} в качестве аппроксимирующей функции [[Интерполяция кубическими сплайнами|кубический сплайн]]:
{{eqno |4}}
{{eqno |4}}
-
::<tex>\int_a^bf(t)dt \approx \sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt,</tex> где
+
::<tex>\int_{t_{i-1}}^{t_i}f(t)dt = \int_{t_{i-1}}^{t_i}\varphi_i(t)dt+R_i,</tex> где
-
<tex>\varphi_i(t) = a_i+b_i(t-t_{i-1})+c_i(t - t_{i-1})^2 + d_i(t - t_{i-1})^3, \; t \in [t_{i-1},\;t_i]. </tex>
+
{{eqno|φ}}
 +
::<tex>\varphi_i(t) = a_i+b_i(t-t_{i-1})+c_i(t - t_{i-1})^2 + d_i(t - t_{i-1})^3, \; t \in [t_{i-1},\;t_i]. </tex>
-
Коэффициенты вычисляются по следующим формулам:
+
Известно, что коэффициенты <tex> a_i, b_i</tex> и <tex>d_i</tex> вычисляются по следующим формулам:
-
{{eqno|5}}
+
{{eqno|5a}}
-
<tex> a_i=f_{i-1}</tex>
+
::<tex> a_i=f_{i-1}</tex>
-
<tex>b_i=\frac{f_i-f_{i-1}}{\tau}\;-\;\frac{\tau(2c_i+c_{i+1})}{3}</tex>
+
{{eqno|5b}}
 +
::<tex>b_i=\frac{f_i-f_{i-1}}{\tau_i}\;-\;\frac{\tau_i(2c_i+c_{i+1})}{3}</tex>
-
<tex>c_{i-1}+4c_i+c_{i+1}=3\left(\frac{f_i-2f_{i-1}+f_{i-2}}{\tau^2}\right); \;\;c_1=0</tex>
+
{{eqno|5d}}
 +
::<tex>d_i=\frac{c_{i+1}-c_i}{3\tau_i}</tex>
-
<tex>d_i=\frac{c_{i+1}-c_i}{3\tau}</tex>
+
А коэффициенты <tex>c_i</tex> являются решением [[СЛАУ]]:
 +
{{eqno|5c}}
 +
::<tex>\tau_{i-1}c_{i-1}+2(\tau_{i-1}+\tau_i)c_i+\tau_i c_{i+1}=3\left(\frac{f_i-f_{i-1}}{\tau_i} \; - \; \frac{f_{i-1}-f_{i-2}}{\tau_{i-1}}\right),\;\; 2 \le i \le N</tex>
 +
Нетрудно видеть, что матрица для решения [[СЛАУ]] {{eqref|5c}} есть трёхдиагональная матрица с диагональным преобладанием. Поэтому коэффициенты <tex>c_i</tex> можно вычислить с помощью метода прогонки.
-
Тогда интеграл {{eqref|4}} запишется как сумма интегралов от сплайнов:
+
Коэффициенты <tex>c_1</tex> и <tex>c_{N+1}</tex> получаются из условий свободных концов сплайна. Обычно требуют нулевую кривизну на концах сплайна и берут <tex>c_1=c_{N+1}=0.</tex>
-
::<tex>J=\int_a^bf(t)dt \approx \sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt=\sum_{i=1}^N\left(a_i\tau+{b_i\over2}\tau^2+{c_i\over3}\tau^3+{d_i\over4}\tau^4\right).</tex>
+
В итоге интеграл {{eqref|1*}} запишется как сумма интегралов от сплайнов:
-
Последняя формула упрощается при подстановке в неё выражений {{eqref|5}} для коэффициентов <tex>a_i,\;b_i, \; d_i:</tex>
+
::<tex>J=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt+R=\sum_{i=1}^N\left(a_i\tau_i+{b_i\over2}\tau_i^2+{c_i\over3}\tau_i^3+{d_i\over4}\tau_i^4\right)+R.</tex>
-
<tex>J \;\approx \;\sum_{i=1}^n\;\frac{f_i+f_{i-1}}{2}\tau\;-\;\sum_{i=1}^n\;\frac{\tau^3(c_{i+1}+c_i)}{12}</tex>
+
Последняя формула упрощается при подстановке в неё выражений {{eqref|5a}}, {{eqref|5b}} и {{eqref|5d}} для коэффициентов <tex>a_i,\;b_i</tex> и <tex>d_i:</tex>
 +
{{eqno|6}}
 +
::<tex>J = \sum_{i=1}^N\;\frac{f_i+f_{i-1}}{2}\tau_i\;-\;\sum_{i=1}^N\;\frac{\tau_i^3(c_{i+1}+c_i)}{12}+R.</tex>
 +
=== Случай с равномерной сеткой ===
 +
Пусть на отрезке <tex>[a,b]</tex> задана равномерная сетка, т.е. <tex>\tau_1=\tau_2=\dots=\tau_N=\tau.</tex> Тогда выражение {{eqref|6}} перепишется в виде:
 +
{{eqno|7}}
 +
::<tex>J\;=\;\frac{\tau}{2}(f_0+f_N)+\tau\sum_{i=1}^{N-1}f_i-\frac{\tau^3}{6}\sum_{i=2}^Nc_i+R.</tex>
 +
 +
Просуммируем уравнения {{eqref|5c}} от i=2 до N. Получим:
 +
 +
{{eqno|8}}
 +
::<tex> 6\sum_{i=2}^Nc_i=(c_2+c_N)+\frac{3}{\tau^2}(f_0-f_1+f_N-f_{N-1}).</tex>
 +
 +
Подставим {{eqref|8}} в {{eqref|7}} и получим окончательное выражение для <tex>J</tex>:
 +
 +
{{eqno|9}}
 +
::<tex> J\;=\;\tau\sum_{i=2}^{N-2}f_i+\frac{\tau}{12}(5f_0+13f_1+13f_{N-1}+5f_N)-\frac{\tau^3}{36}(c_2+c_N)+R</tex>
 +
 +
Несмотря на то, что <tex>c_2</tex> и <tex>c_N</tex> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
== Анализ метода и ошибок ==
== Анализ метода и ошибок ==
-
<!--у Буслова про поправочный член в трапециях-->
+
 
 +
Анализ формулы {{eqref|6}} показывает, что первый член в правой части совпадает с [[Методы прямоугольников и трапеций|формулой трапеций]]. Следовательно, второй член характеризует поправку к методу трапеций, которую дает использование сплайнов.
 +
 
 +
Как следует из формулы {{eqref|φ}}, коэффициенты <tex>c_i</tex> выражаются через вторые производные <tex>\varphi_i''(t):</tex>
 +
 
 +
::<tex>c_i=\frac{1}{2}\varphi_i''(t_{i-1})\approx\frac{1}{2}f_{i-1}''.</tex>
 +
 
 +
Это позволяет оценить второй член правой части формулы {{eqref|6}}:
 +
 
 +
::<tex>\frac{\tau_i^3}{12}(c_{i-1}+c_i)\approx\frac{\tau_i^3}{12}f_*'',</tex>
 +
 
 +
где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
 +
 
 +
===Ошибки округления===
 +
 
 +
Следует помнить, что при реализации метода помимо теоретической априорной погрешности из-за суммы большого числа чисел возникает также погрешность округления.
 +
 
== Числовой пример ==
== Числовой пример ==
 +
 +
[[Изображение:Splineint.jpg|thumb|200px| рис.1. График функции]]
 +
Рассмотрим функцию <tex>f(t)=sin(t).</tex> Вычислим с помощью сплайн-квадратур приближенное значение интеграла
 +
 +
::<tex>\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2</tex>
 +
 +
и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:
 +
 +
{| class="wikitable"
 +
|[[Изображение:Splineint100.jpg|thumb|400px| рис.2. Число разбиений: 100]]
 +
|[[Изображение:Splineint1kk.jpg|thumb|400px| рис.3. Число разбиений: 1 000 000]]
 +
|}
 +
 +
При этом красным цветом показан результат рассчетов по формуле {{eqref|9}}, а зеленым - по формуле {{eqref|6}}.
 +
 +
Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.
 +
== Рекомендации программисту ==
== Рекомендации программисту ==
-
<!-- У Богдана интересная хрень про заданную точность -->
+
 
 +
=== Пример программы ===
 +
 
 +
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
 +
[[Media:Splineint.zip|Splineint.zip [0.7Кб]]]
 +
 
 +
'''Некоторые комментарии по работе с программой:'''
 +
 
 +
В 6-й строке <code>const int N=1000;</code> N - число отрезков <tex>\tau_i.</tex>
 +
 
 +
В 8-й строке <code>const double a=0,b=pi;</code> <tex>a</tex> и <tex>b</tex> - пределы интегрирования.
 +
 
 +
В 18-й строке <code>f[i]=sin(x);</code> f[i] - интегрируемая функция.
 +
 
== Заключение ==
== Заключение ==
-
<!-- В Мудрове есть че-то на 158 стр -->
+
 
 +
Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако [[алгоритм]] интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <tex>c_i.</tex> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны уже используются, в частности, для сглаживания зависимостей, обработки эксперимениальных данных и т.п.
 +
 
 +
Также одним из основных преимуществ метода является возможность работать с произвольным разбиением орезка.
 +
 
== Ссылки ==
== Ссылки ==
* [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
* [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
 +
* [[Интерполяция кубическими сплайнами]]
 +
== Список литературы ==
== Список литературы ==
* http://www.intuit.ru/department/calculate/calcmathbase/7/1.html
* http://www.intuit.ru/department/calculate/calcmathbase/7/1.html
* http://mathalgo.blogspot.com/2007/11/blog-post.html
* http://mathalgo.blogspot.com/2007/11/blog-post.html
-
* ''Бабенко К.И'' Основы численного анализа М.: Наука, 1986.
+
* http://myhomepage.h17.ru/Lect06/lect06.htm#02
-
* Мудров
+
* А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.
-
* Буслов
+
-
<!--не забыть про ссылки и норм. лит-ру -->
 
-
{{stub}}
 
[[Категория:Численное интегрирование]]
[[Категория:Численное интегрирование]]

Текущая версия

Содержание

Введение

Ставится задача вычислить интеграл вида

(1)
J=\int_a^bf(t)dt,

где a и b - нижний и верхний пределы интегрирования; f(t) - непрерывная функция на отрезке [a,b].

Пусть на отрезке интегрирования дана сетка, т.е. имеется совокупность узлов \left\{{t_i}\right\}_{i = 0}^{N},\; i\in \left[{a, b}\right]. Тогда интервал [a,b] разобьется на участки \tau_i = t_i-t_{i-1},\; i=1,\dots,N.. Пусть также известны значения функции в узлах этой сетки, т.е. задана таблица f_i = \left\{{f(t_i)}\right\}_{i = 0}^{N}. Представим интеграл (1) в виде суммы интегралов по частичным отрезкам:

(2)
\int_a^bf(t)dt=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}f(t)dt.

Сущность большинства методов вычисления определённых интегралов состоит в замене подынтегральной функции f(t) на отрезке [t_{i-1},\;t_i] аппроксимирующей функцией \varphi(t), для которой можно легко записать первообразную в элементарных функциях, т.е.

(3)
\int_{t_{i-1}}^{t_i}f(t)dt=\int_{t_{i-1}}^{t_i}\varphi(t)dt+R_i=S_i+R_i,
(1*)
J=\sum_{i=1}^N \left( \int_{t_{i-1}}^{t_i}\varphi(t)dt+R_i \right)=S+R,

где S_i - приближённое значение интеграла на i-м отрезке, S=\sum_{i=1}^NS_i, а R_i - погрешность вычисления интеграла, R=\sum_{i=1}^NR_i,. Лучше всего изучена замена f(t) алгебраическим многочленом.

Изложение метода

Общий случай

Возьмём в (3) в качестве аппроксимирующей функции кубический сплайн:

(4)
\int_{t_{i-1}}^{t_i}f(t)dt = \int_{t_{i-1}}^{t_i}\varphi_i(t)dt+R_i, где
(φ)
\varphi_i(t) = a_i+b_i(t-t_{i-1})+c_i(t - t_{i-1})^2 + d_i(t - t_{i-1})^3, \; t \in [t_{i-1},\;t_i].

Известно, что коэффициенты  a_i, b_i и d_i вычисляются по следующим формулам:

(5a)
 a_i=f_{i-1}
(5b)
b_i=\frac{f_i-f_{i-1}}{\tau_i}\;-\;\frac{\tau_i(2c_i+c_{i+1})}{3}
(5d)
d_i=\frac{c_{i+1}-c_i}{3\tau_i}

А коэффициенты c_i являются решением СЛАУ:

(5c)
\tau_{i-1}c_{i-1}+2(\tau_{i-1}+\tau_i)c_i+\tau_i c_{i+1}=3\left(\frac{f_i-f_{i-1}}{\tau_i} \; - \; \frac{f_{i-1}-f_{i-2}}{\tau_{i-1}}\right),\;\; 2 \le i \le N

Нетрудно видеть, что матрица для решения СЛАУ (5c) есть трёхдиагональная матрица с диагональным преобладанием. Поэтому коэффициенты c_i можно вычислить с помощью метода прогонки.

Коэффициенты c_1 и c_{N+1} получаются из условий свободных концов сплайна. Обычно требуют нулевую кривизну на концах сплайна и берут c_1=c_{N+1}=0.

В итоге интеграл (1*) запишется как сумма интегралов от сплайнов:

J=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt+R=\sum_{i=1}^N\left(a_i\tau_i+{b_i\over2}\tau_i^2+{c_i\over3}\tau_i^3+{d_i\over4}\tau_i^4\right)+R.

Последняя формула упрощается при подстановке в неё выражений (5a), (5b) и (5d) для коэффициентов a_i,\;b_i и d_i:

(6)
J = \sum_{i=1}^N\;\frac{f_i+f_{i-1}}{2}\tau_i\;-\;\sum_{i=1}^N\;\frac{\tau_i^3(c_{i+1}+c_i)}{12}+R.

Случай с равномерной сеткой

Пусть на отрезке [a,b] задана равномерная сетка, т.е. \tau_1=\tau_2=\dots=\tau_N=\tau. Тогда выражение (6) перепишется в виде:

(7)
J\;=\;\frac{\tau}{2}(f_0+f_N)+\tau\sum_{i=1}^{N-1}f_i-\frac{\tau^3}{6}\sum_{i=2}^Nc_i+R.

Просуммируем уравнения (5c) от i=2 до N. Получим:

(8)
 6\sum_{i=2}^Nc_i=(c_2+c_N)+\frac{3}{\tau^2}(f_0-f_1+f_N-f_{N-1}).

Подставим (8) в (7) и получим окончательное выражение для J:

(9)
 J\;=\;\tau\sum_{i=2}^{N-2}f_i+\frac{\tau}{12}(5f_0+13f_1+13f_{N-1}+5f_N)-\frac{\tau^3}{36}(c_2+c_N)+R

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

Анализ метода и ошибок

Анализ формулы (6) показывает, что первый член в правой части совпадает с формулой трапеций. Следовательно, второй член характеризует поправку к методу трапеций, которую дает использование сплайнов.

Как следует из формулы (φ), коэффициенты c_i выражаются через вторые производные \varphi_i''(t):

c_i=\frac{1}{2}\varphi_i''(t_{i-1})\approx\frac{1}{2}f_{i-1}''.

Это позволяет оценить второй член правой части формулы (6):

\frac{\tau_i^3}{12}(c_{i-1}+c_i)\approx\frac{\tau_i^3}{12}f_*'',

где f_*'' - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.

Ошибки округления

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

Числовой пример

рис.1. График функции
рис.1. График функции

Рассмотрим функцию f(t)=sin(t). Вычислим с помощью сплайн-квадратур приближенное значение интеграла

\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2

и исследуем поведение погрешности. Результаты работы программы проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:

рис.2. Число разбиений: 100
рис.2. Число разбиений: 100
рис.3. Число разбиений: 1 000 000
рис.3. Число разбиений: 1 000 000

При этом красным цветом показан результат рассчетов по формуле (9), а зеленым - по формуле (6).

Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.

Рекомендации программисту

Пример программы

Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур: Splineint.zip [0.7Кб]

Некоторые комментарии по работе с программой:

В 6-й строке const int N=1000; N - число отрезков \tau_i.

В 8-й строке const double a=0,b=pi; a и b - пределы интегрирования.

В 18-й строке f[i]=sin(x); f[i] - интегрируемая функция.

Заключение

Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако алгоритм интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов c_i. Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны уже используются, в частности, для сглаживания зависимостей, обработки эксперимениальных данных и т.п.

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

Ссылки

Список литературы

Личные инструменты