Участник:Александр Двойнев/песочница

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

(Различия между версиями)
Перейти к: навигация, поиск
Строка 60: Строка 60:
{{eqno|7}}
{{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>
+
::<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. Получим:
Просуммируем уравнения {{eqref|5c}} от i=2 до N. Получим:
Строка 69: Строка 69:
Подставим {{eqref|8}} в {{eqref|7}} и получим окончательное выражение для <tex>J</tex>:
Подставим {{eqref|8}} в {{eqref|7}} и получим окончательное выражение для <tex>J</tex>:
-
::<tex> J\;\approx\;\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)</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> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
Несмотря на то, что <tex>c_2</tex> и <tex>c_N</tex> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
Строка 89: Строка 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}}.
== Рекомендации программисту ==
== Рекомендации программисту ==

Версия 14:01, 27 ноября 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_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).

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

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

Ниже приведен пример программы на языке 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] - интегрируемая функция.

Заключение

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

Ссылки

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

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%A3%D1%87%D0%B0%D1%81%D1%82%D0%BD%D0%B8%D0%BA:%D0%90%D0%BB%D0%B5%D0%BA%D1%81%D0%B0%D0%BD%D0%B4%D1%80_%D0%94%D0%B2%D0%BE%D0%B9%D0%BD%D0%B5%D0%B2/%D0%BF%D0%B5%D1%81%D0%BE%D1%87%D0%BD%D0%B8%D1%86%D0%B0»
Личные инструменты