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

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

(Различия между версиями)
Перейти к: навигация, поиск
(Добавил основную часть теории)
Строка 6: Строка 6:
где <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>
Строка 20: Строка 20:
== Изложение метода ==
== Изложение метода ==
-
 
+
Возьмем в {{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_a^bf(t)dt \approx \sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt,</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>
Коэффициенты вычисляются по следующим формулам:
Коэффициенты вычисляются по следующим формулам:
-
{{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|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, \;\;c_1=c_{N+1}=0</tex>
-
<tex>d_i=\frac{c_{i+1}-c_i}{3\tau}</tex>
+
{{eqno|5d}}
 +
::<tex>d_i=\frac{c_{i+1}-c_i}{3\tau_i}</tex>
Тогда интеграл {{eqref|4}} запишется как сумма интегралов от сплайнов:
Тогда интеграл {{eqref|4}} запишется как сумма интегралов от сплайнов:
-
::<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>
+
::<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_i+{b_i\over2}\tau_i^2+{c_i\over3}\tau_i^3+{d_i\over4}\tau_i^4\right).</tex>
-
Последняя формула упрощается при подстановке в неё выражений {{eqref|5}} для коэффициентов <tex>a_i,\;b_i, \; d_i:</tex>
+
Последняя формула упрощается при подстановке в неё выражений {{eqref|5a}}, {{eqref|5b}} и {{eqref|5d}} для коэффициентов <tex>a_i,\;b_i</tex> и <tex>d_i:</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>
+
{{eqno|6}}
 +
::<tex>J \;\approx \;\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}.</tex>
 +
Нетрудно видеть, что матриа для решения СЛАУ {{eqref|5c}} есть трехдиагональная матрица с диагональным преобладанием. Поэтому коэффициенты <tex>c_i</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> - вторая производная в некоторой нутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
-
== Анализ метода и ошибок ==
 
-
<!--у Буслова про поправочный член в трапециях-->
 
== Числовой пример ==
== Числовой пример ==
 +
 +
== Рекомендации программисту ==
== Рекомендации программисту ==
 +
===Пример программы===
 +
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
 +
[[Media:Splineint.zip|Splineint.zip [0.7Кб]]]
 +
 +
'''Некоторые комментарии по работе с программой:'''
 +
 +
В 5-й строке <code>const int N=100;</code> N - число отрезков <tex>\tau_i.</tex>
 +
 +
В 7-й строке <code>const double a=1,b=6;</code> <tex>a</tex> и <tex>b</tex> - пределы интегрирования.
 +
 +
В 49-й строке <code>f[i]=0.6*x*x*x-3*x*x+3*x;</code> f[i] - интегрируемая функция.
 +
 +
 +
=== Случай с равномерной сеткой ===
 +
Пусть на отрезке <tex>[a,b]</tex> задана равномерная сетка, т.е. <tex>\tau_1=\tau_2=\dots=\tau_N=\tau.</tex> Тогда выражение {{eqref|6}} перепишется в виде:
 +
 +
{{eqno|7}}
 +
::<tex>J\;\approx\;\frac{\tau}{2}(f_0+f_N)+\tau\sum_{i=1}^{N-1}f_i-\frac{\tau^3}{6}\sum_{i=2}^Nc_i.</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>:
 +
 +
::<tex> J\;\approx\;\tau\sum_{i=2}^{N-2}f_i+\frac{\tau}{2}(-5f_0+8f_1+8f_{N-1}-5f_N)-\tau^3(c_2+c_N)</tex>
 +
 +
Несмотря на то, что <tex>c_2</tex> и <tex>c_N</tex> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
 +
<!-- У Богдана интересная хрень про заданную точность -->
<!-- У Богдана интересная хрень про заданную точность -->
== Заключение ==
== Заключение ==
-
<!-- В Мудрове есть че-то на 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}}
{{stub}}
[[Категория:Численное интегрирование]]
[[Категория:Численное интегрирование]]

Версия 10:02, 19 октября 2008

Содержание

Введение

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

(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=S+R,

где S - приближеное значение интеграла; R - погрешность вычисления интеграла. Лучше всего изучена замена f(t) алгебраическим многочленом.

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

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

(4)
\int_a^bf(t)dt \approx \sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt, где
(φ)
\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].

Коэффициенты вычисляются по следующим формулам:

(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}
(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, \;\;c_1=c_{N+1}=0
(5d)
d_i=\frac{c_{i+1}-c_i}{3\tau_i}

Тогда интеграл (4) запишется как сумма интегралов от сплайнов:

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_i+{b_i\over2}\tau_i^2+{c_i\over3}\tau_i^3+{d_i\over4}\tau_i^4\right).

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

(6)
J \;\approx \;\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}.

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

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

Анализ формулы (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_*'' - вторая производная в некоторой нутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.

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

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

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

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

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

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

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

В 49-й строке f[i]=0.6*x*x*x-3*x*x+3*x; f[i] - интегрируемая функция.


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

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

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

Просуммируем уравнения (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:

J\;\approx\;\tau\sum_{i=2}^{N-2}f_i+\frac{\tau}{2}(-5f_0+8f_1+8f_{N-1}-5f_N)-\tau^3(c_2+c_N)

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

Заключение

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

Ссылки

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

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%9F%D1%80%D0%B8%D0%BC%D0%B5%D0%BD%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D0%BF%D0%BB%D0%B0%D0%B9%D0%BD%D0%BE%D0%B2_%D0%B4%D0%BB%D1%8F_%D1%87%D0%B8%D1%81%D0%BB%D0%B5%D0%BD%D0%BD%D0%BE%D0%B3%D0%BE_%D0%B8%D0%BD%D1%82%D0%B5%D0%B3%D1%80%D0%B8%D1%80%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D1%8F»
Личные инструменты