Применение сплайнов для численного интегрирования
Материал из MachineLearning.
(орфография, категория) |
(категория) |
||
(4 промежуточные версии не показаны) | |||
Строка 5: | Строка 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 | + | <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}} в виде суммы интегралов по частичным отрезкам: | ||
Строка 15: | Строка 15: | ||
Сущность большинства методов вычисления определённых интегралов состоит в замене подынтегральной функции <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+ | + | ::<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> |
- | где | + | {{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>\ | + | ::<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> | ||
Строка 32: | Строка 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|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> | ||
- | + | А коэффициенты <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= | + | ::<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 | + | ::<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> |
- | + | === Случай с равномерной сеткой === | |
+ | |||
+ | Пусть на отрезке <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> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых. | ||
== Анализ метода и ошибок == | == Анализ метода и ошибок == | ||
Строка 63: | Строка 87: | ||
где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций. | где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций. | ||
+ | |||
+ | ===Ошибки округления=== | ||
+ | |||
+ | Следует помнить, что при реализации метода помимо теоретической априорной погрешности из-за суммы большого числа чисел возникает также погрешность округления. | ||
== Числовой пример == | == Числовой пример == | ||
- | [[Изображение:Splineint.jpg|thumb|200px| График функции]] | + | [[Изображение:Splineint.jpg|thumb|200px| рис.1. График функции]] |
- | Рассмотрим функцию <tex>f(t)= | + | Рассмотрим функцию <tex>f(t)=sin(t).</tex> Вычислим с помощью сплайн-квадратур приближенное значение интеграла |
- | ::<tex>\int_0^ | + | ::<tex>\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2</tex> |
- | и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] | + | и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений: |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | Как видно из | + | {| class="wikitable" |
+ | |[[Изображение:Splineint100.jpg|thumb|400px| рис.2. Число разбиений: 100]] | ||
+ | |[[Изображение:Splineint1kk.jpg|thumb|400px| рис.3. Число разбиений: 1 000 000]] | ||
+ | |} | ||
+ | |||
+ | При этом красным цветом показан результат рассчетов по формуле {{eqref|9}}, а зеленым - по формуле {{eqref|6}}. | ||
+ | |||
+ | Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти. | ||
== Рекомендации программисту == | == Рекомендации программисту == | ||
- | ===Пример программы=== | + | === Пример программы === |
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур: | Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур: | ||
Строка 127: | Строка 119: | ||
'''Некоторые комментарии по работе с программой:''' | '''Некоторые комментарии по работе с программой:''' | ||
- | В | + | В 6-й строке <code>const int N=1000;</code> N - число отрезков <tex>\tau_i.</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> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны уже используются, в частности, для сглаживания зависимостей, обработки эксперимениальных данных и т.п. | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | Также одним из основных преимуществ метода является возможность работать с произвольным разбиением орезка. | |
== Ссылки == | == Ссылки == | ||
Строка 166: | Строка 142: | ||
[[Категория:Численное интегрирование]] | [[Категория:Численное интегрирование]] | ||
- |
Текущая версия
Содержание |
Введение
Ставится задача вычислить интеграл вида
где и - нижний и верхний пределы интегрирования; - непрерывная функция на отрезке .
Пусть на отрезке интегрирования дана сетка, т.е. имеется совокупность узлов Тогда интервал разобьется на участки . Пусть также известны значения функции в узлах этой сетки, т.е. задана таблица Представим интеграл (1) в виде суммы интегралов по частичным отрезкам:
Сущность большинства методов вычисления определённых интегралов состоит в замене подынтегральной функции на отрезке аппроксимирующей функцией , для которой можно легко записать первообразную в элементарных функциях, т.е.
где - приближённое значение интеграла на i-м отрезке, а - погрешность вычисления интеграла, . Лучше всего изучена замена алгебраическим многочленом.
Изложение метода
Общий случай
Возьмём в (3) в качестве аппроксимирующей функции кубический сплайн:
- где
Известно, что коэффициенты и вычисляются по следующим формулам:
А коэффициенты являются решением СЛАУ:
Нетрудно видеть, что матрица для решения СЛАУ (5c) есть трёхдиагональная матрица с диагональным преобладанием. Поэтому коэффициенты можно вычислить с помощью метода прогонки.
Коэффициенты и получаются из условий свободных концов сплайна. Обычно требуют нулевую кривизну на концах сплайна и берут
В итоге интеграл (1*) запишется как сумма интегралов от сплайнов:
Последняя формула упрощается при подстановке в неё выражений (5a), (5b) и (5d) для коэффициентов и
Случай с равномерной сеткой
Пусть на отрезке задана равномерная сетка, т.е. Тогда выражение (6) перепишется в виде:
Просуммируем уравнения (5c) от i=2 до N. Получим:
Подставим (8) в (7) и получим окончательное выражение для :
Несмотря на то, что и все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
Анализ метода и ошибок
Анализ формулы (6) показывает, что первый член в правой части совпадает с формулой трапеций. Следовательно, второй член характеризует поправку к методу трапеций, которую дает использование сплайнов.
Как следует из формулы (φ), коэффициенты выражаются через вторые производные
Это позволяет оценить второй член правой части формулы (6):
где - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
Ошибки округления
Следует помнить, что при реализации метода помимо теоретической априорной погрешности из-за суммы большого числа чисел возникает также погрешность округления.
Числовой пример
Рассмотрим функцию Вычислим с помощью сплайн-квадратур приближенное значение интеграла
и исследуем поведение погрешности. Результаты работы программы проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:
При этом красным цветом показан результат рассчетов по формуле (9), а зеленым - по формуле (6).
Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.
Рекомендации программисту
Пример программы
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур: Splineint.zip [0.7Кб]
Некоторые комментарии по работе с программой:
В 6-й строке const int N=1000;
N - число отрезков
В 8-й строке const double a=0,b=pi;
и - пределы интегрирования.
В 18-й строке f[i]=sin(x);
f[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.