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