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

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

(Различия между версиями)
Перейти к: навигация, поиск
(категория)
 
(5 промежуточных версий не показаны.)
Строка 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>c_1</tex> и <tex>c_{N+1}</tex> получаются из условий свободных концов сплайна. Обычно требуют нулевую кривизну на концах сплайна и берут <tex>c_1=c_{N+1}=0.</tex>
 +
 
 +
В итоге интеграл {{eqref|1*}} запишется как сумма интегралов от сплайнов:
-
::<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>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> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
== Анализ метода и ошибок ==
== Анализ метода и ошибок ==
Строка 64: Строка 87:
где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
 +
 +
===Ошибки округления===
 +
 +
Следует помнить, что при реализации метода помимо теоретической априорной погрешности из-за суммы большого числа чисел возникает также погрешность округления.
== Числовой пример ==
== Числовой пример ==
-
[[Изображение:Splineint.jpg|thumb|200px| График функции]]
+
[[Изображение:Splineint.jpg|thumb|200px| рис.1. График функции]]
-
Рассмотрим функцию <tex>f(t)=0.6t^3-3t^2+3t.</tex> Вычислим с помощью сплайн-квадратур приближенное значение интеграла
+
Рассмотрим функцию <tex>f(t)=sin(t).</tex> Вычислим с помощью сплайн-квадратур приближенное значение интеграла
 +
 
 +
::<tex>\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2</tex>
 +
 
 +
и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:
-
::<tex>\int_0^4f(t)dt=\int_0^4(0.6t^3-3t^2+3t)dt=-1,6</tex>
+
{| class="wikitable"
 +
|[[Изображение:Splineint100.jpg|thumb|400px| рис.2. Число разбиений: 100]]
 +
|[[Изображение:Splineint1kk.jpg|thumb|400px| рис.3. Число разбиений: 1 000 000]]
 +
|}
-
и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] приведены в следующей таблице:
+
При этом красным цветом показан результат рассчетов по формуле {{eqref|9}}, а зеленым - по формуле {{eqref|6}}.
-
+
-
::{| class="standard"
+
-
!N
+
-
!J
+
-
!ε
+
-
|-
+
-
!5
+
-
| -8.88236
+
-
|7.28236
+
-
|-
+
-
!10
+
-
| -3.61479
+
-
|2.01479
+
-
|-
+
-
!20
+
-
| -2.13136
+
-
|0.53136
+
-
|-
+
-
!50
+
-
| -1.68776
+
-
|0.087758
+
-
|-
+
-
!100
+
-
| -1.62217
+
-
|0.022169
+
-
|-
+
-
!200
+
-
| -1.60557
+
-
|0.00557
+
-
|-
+
-
!500
+
-
| -1.60089
+
-
|0.00089
+
-
|-
+
-
!1000
+
-
| -1.60022
+
-
|0.00022
+
-
|-
+
-
!2000
+
-
| -1.60006
+
-
|0.00005
+
-
|}
+
-
Здесь N - число отрезков, на которые разбивается интервал [0,4], J - приблизительное значение интеграла, ε - ошибка.
+
-
Как видно из таблицы, при малых N (особенно при N=5) ошибка невероятно велика. Однако с ростом N ошибка стремительно убывает, и приблизительное значение интеграла сходится к правильному значению.
+
Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.
== Рекомендации программисту ==
== Рекомендации программисту ==
-
===Пример программы===
+
=== Пример программы ===
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
Строка 128: Строка 119:
'''Некоторые комментарии по работе с программой:'''
'''Некоторые комментарии по работе с программой:'''
-
В 5-й строке <code>const int N=100;</code> N - число отрезков <tex>\tau_i.</tex>
+
В 6-й строке <code>const int N=1000;</code> N - число отрезков <tex>\tau_i.</tex>
-
В 7-й строке <code>const double a=1,b=6;</code> <tex>a</tex> и <tex>b</tex> - пределы интегрирования.
+
В 8-й строке <code>const double a=0,b=pi;</code> <tex>a</tex> и <tex>b</tex> - пределы интегрирования.
-
В 49-й строке <code>f[i]=0.6*x*x*x-3*x*x+3*x;</code> f[i] - интегрируемая функция.
+
В 18-й строке <code>f[i]=sin(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> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
+
== Заключение ==
== Заключение ==
-
Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако алгоритм интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <tex>c_i.</tex> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны применяются для сглаживания зависимостей, обработки эксперимениальных данных и т.п.
+
Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако [[алгоритм]] интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <tex>c_i.</tex> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны уже используются, в частности, для сглаживания зависимостей, обработки эксперимениальных данных и т.п.
 +
 
 +
Также одним из основных преимуществ метода является возможность работать с произвольным разбиением орезка.
== Ссылки ==
== Ссылки ==
Строка 165: Строка 140:
* http://myhomepage.h17.ru/Lect06/lect06.htm#02
* http://myhomepage.h17.ru/Lect06/lect06.htm#02
* А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.
* А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.
 +
[[Категория:Численное интегрирование]]
[[Категория:Численное интегрирование]]

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

Содержание

Введение

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

(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»
Личные инструменты