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

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

(Различия между версиями)

Перейти к: навигация, поиск

Версия 14:54, 30 ноября 2008

Содержание

1 Введение
2 Изложение метода
- 2.1 Общий случай
- 2.2 Случай с равномерной сеткой
3 Анализ метода и ошибок
4 Числовой пример
5 Рекомендации программисту
- 5.1 Пример программы
- 5.2 Ошибки округления
6 Заключение
7 Ссылки
8 Список литературы

Введение

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

(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. График функции

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

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

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

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

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

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

Заключение

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

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

Ссылки

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

http://www.intuit.ru/department/calculate/calcmathbase/7/1.html
http://mathalgo.blogspot.com/2007/11/blog-post.html
http://myhomepage.h17.ru/Lect06/lect06.htm#02
А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.

Источник — «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»

@@ Строка 5: / Строка 5: @@
 где <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>
 Представим интеграл {{eqref|1}} в виде суммы интегралов по частичным отрезкам:
@@ Строка 15: / Строка 15: @@
 Сущность большинства методов вычисления определённых интегралов состоит в замене подынтегральной функции <tex>f(t)</tex> на отрезке <tex>[t_{i-1},\;t_i]</tex> аппроксимирующей функцией <tex>\varphi(t)</tex>, для которой можно легко записать первообразную в элементарных функциях, т.е.
 {{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}} в качестве аппроксимирующей функции [[Интерполяция кубическими сплайнами|кубический сплайн]]:
 {{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|φ}}
 ::<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}}
 ::<tex> a_i=f_{i-1}</tex>
@@ Строка 32: / Строка 35: @@
 {{eqno|5b}}
 ::<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}}
 ::<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>
 {{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> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
 == Анализ метода и ошибок ==
@@ Строка 66: / Строка 90: @@
 == Числовой пример ==
-[[Изображение: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^4f(t)dt=\int_0^4(0.6t^3-3t^2+3t)dt=-1,6</tex>
+::<tex>\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2</tex>
-и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] приведены в следующей таблице:
+и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:
-::{| 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 ошибка стремительно убывает, и приблизительное значение интеграла сходится к правильному значению.
+{| class="wikitable"
+|[[Изображение:Splineint100.jpg|thumb|400px| рис.2. Число разбиений: 100]]
+|[[Изображение:Splineint1kk.jpg|thumb|400px| рис.3. Число разбиений: 1 000 000]]
+|}
+При этом красным цветом показан результат рассчетов по формуле {{eqref|9}}, а зеленым - по формуле {{eqref|6}}.
+Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.
 == Рекомендации программисту ==
-===Пример программы===
+=== Пример программы ===
 Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
@@ Строка 127: / Строка 115: @@
 '''Некоторые комментарии по работе с программой:'''
-В 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}{12}(5f_0+7f_1+5f_{N-1}+7f_N)-\frac{\tau^3}{36}(c_2+c_N)</tex>
-Несмотря на то, что <tex>c_2</tex> и <tex>c_N</tex> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
 == Заключение ==
 Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако [[алгоритм]] интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <tex>c_i.</tex> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны применяются для сглаживания зависимостей, обработки эксперимениальных данных и т.п.
+Также одним из основных преимуществ метода является возможность работать с произвольным разбиением орезка.
 == Ссылки ==
@@ Строка 164: / Строка 140: @@
 * http://myhomepage.h17.ru/Lect06/lect06.htm#02
 * А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.
-[[Категория:Численное интегрирование]]
-[[Категория:Учебные задачи]]

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

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

Версия 14:54, 30 ноября 2008

Содержание

Введение

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

Общий случай

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

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

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

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

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

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

Заключение

Ссылки

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

Просмотры

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

Навигация

Поиск

Инструменты