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

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

(Различия между версиями)
Перейти к: навигация, поиск
(категория)
 
(6 промежуточных версий не показаны.)
Строка 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},\; 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>\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> где
{{eqno|φ}}
{{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>\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|5a}}
{{eqno|5a}}
::<tex> a_i=f_{i-1}</tex>
::<tex> a_i=f_{i-1}</tex>
Строка 33: Строка 35:
{{eqno|5b}}
{{eqno|5b}}
::<tex>b_i=\frac{f_i-f_{i-1}}{\tau_i}\;-\;\frac{\tau_i(2c_i+c_{i+1})}{3}</tex>
::<tex>b_i=\frac{f_i-f_{i-1}}{\tau_i}\;-\;\frac{\tau_i(2c_i+c_{i+1})}{3}</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>
 
{{eqno|5d}}
{{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_i}</tex>
-
Тогда интеграл {{eqref|4}} запишется как сумма интегралов от сплайнов:
+
А коэффициенты <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> можно вычислить с помощью метода прогонки.
-
::<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>
+
Коэффициенты <tex>c_1</tex> и <tex>c_{N+1}</tex> получаются из условий свободных концов сплайна. Обычно требуют нулевую кривизну на концах сплайна и берут <tex>c_1=c_{N+1}=0.</tex>
 +
 
 +
В итоге интеграл {{eqref|1*}} запишется как сумма интегралов от сплайнов:
 +
 
 +
::<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>
Последняя формула упрощается при подстановке в неё выражений {{eqref|5a}}, {{eqref|5b}} и {{eqref|5d}} для коэффициентов <tex>a_i,\;b_i</tex> и <tex>d_i:</tex>
Последняя формула упрощается при подстановке в неё выражений {{eqref|5a}}, {{eqref|5b}} и {{eqref|5d}} для коэффициентов <tex>a_i,\;b_i</tex> и <tex>d_i:</tex>
{{eqno|6}}
{{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>
+
::<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>
-
Нетрудно видеть, что матриа для решения СЛАУ {{eqref|5c}} есть трехдиагональная матрица с диагональным преобладанием. Поэтому коэффициенты <tex>c_i</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|6}} показывает, что первый член в правой части совпадает с [[Методы прямоугольников и трапеций|формулой трапеций]]. Следовательно, второй член характеризует поправку к методу трапеций, которую дает использование сплайнов.
Строка 62: Строка 86:
::<tex>\frac{\tau_i^3}{12}(c_{i-1}+c_i)\approx\frac{\tau_i^3}{12}f_*'',</tex>
::<tex>\frac{\tau_i^3}{12}(c_{i-1}+c_i)\approx\frac{\tau_i^3}{12}f_*'',</tex>
-
где <tex>f_*''</tex> - вторая производная в некоторой нутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
+
где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
-
== Числовой пример ==
+
===Ошибки округления===
 +
Следует помнить, что при реализации метода помимо теоретической априорной погрешности из-за суммы большого числа чисел возникает также погрешность округления.
-
== Рекомендации программисту ==
+
== Числовой пример ==
-
===Пример программы===
+
-
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
+
-
[[Media:Splineint.zip|Splineint.zip [0.7Кб]]]
+
-
'''Некоторые комментарии по работе с программой:'''
+
[[Изображение:Splineint.jpg|thumb|200px| рис.1. График функции]]
 +
Рассмотрим функцию <tex>f(t)=sin(t).</tex> Вычислим с помощью сплайн-квадратур приближенное значение интеграла
-
В 5-й строке <code>const int N=100;</code> N - число отрезков <tex>\tau_i.</tex>
+
::<tex>\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2</tex>
-
В 7-й строке <code>const double a=1,b=6;</code> <tex>a</tex> и <tex>b</tex> - пределы интегрирования.
+
и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:
-
В 49-й строке <code>f[i]=0.6*x*x*x-3*x*x+3*x;</code> f[i] - интегрируемая функция.
+
{| class="wikitable"
 +
|[[Изображение:Splineint100.jpg|thumb|400px| рис.2. Число разбиений: 100]]
 +
|[[Изображение:Splineint1kk.jpg|thumb|400px| рис.3. Число разбиений: 1 000 000]]
 +
|}
 +
При этом красным цветом показан результат рассчетов по формуле {{eqref|9}}, а зеленым - по формуле {{eqref|6}}.
-
=== Случай с равномерной сеткой ===
+
Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.
-
Пусть на отрезке <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}}
+
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
-
::<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>
+
[[Media:Splineint.zip|Splineint.zip [0.7Кб]]]
-
Подставим {{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>
+
В 6-й строке <code>const int N=1000;</code> N - число отрезков <tex>\tau_i.</tex>
-
Несмотря на то, что <tex>c_2</tex> и <tex>c_N</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] - интегрируемая функция.
-
<!-- У Богдана интересная хрень про заданную точность -->
 
== Заключение ==
== Заключение ==
-
Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако алгоритм интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <tex>c_i.</tex> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны применяются для сглаживания зависимостей, обработки эксперимениальных данных и т.п.
+
 
 +
Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако [[алгоритм]] интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <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
Строка 110: Строка 141:
* А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.
* А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 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. Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны уже используются, в частности, для сглаживания зависимостей, обработки эксперимениальных данных и т.п.

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

Ссылки

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

Источник — «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»
Личные инструменты