Экстраполяция Ричардсона, оценки по Рунге и Эйткену, вычисление интегралов с заданной точностью
Материал из MachineLearning.
(Новая: == Введение == === Постановка математической задачи === == Изложение метода == == Анализ метода == == Числово...) |
(→Заключение) |
||
(17 промежуточных версий не показаны.) | |||
Строка 2: | Строка 2: | ||
=== Постановка математической задачи === | === Постановка математической задачи === | ||
+ | |||
+ | Задача численного интегрирования состоит в приближенном нахождении значения интеграла | ||
+ | {{ eqno | 1}} | ||
+ | :<tex>I = \int\limits_a^b f(x)\,dx,</tex> | ||
+ | |||
+ | где | ||
+ | <tex> | ||
+ | f(x) | ||
+ | </tex> - заданная и интегрируемая на <tex> [a, b] </tex> функция. В качестве приближенного значения рассматривается число | ||
+ | |||
+ | {{ eqno | 2}} | ||
+ | :<tex>I_n=\sum_{i=0}^n c_k f(x_k),</tex> | ||
+ | |||
+ | где <tex>c_k</tex> - числовые коэффициенты и <tex>x_k</tex> - точки отрезка <tex>[a,b]</tex>, <tex> k = 0, 1, \ldots, n </tex>. | ||
+ | Приближенное равенство | ||
+ | |||
+ | :<tex>\int\limits_a^b f(x)\,dx=\sum_{k=0}^n c_k f(x_k)</tex> | ||
+ | |||
+ | называется <i>квадратурной формулой</i>, а сумма вида {{eqref|2}} - <i>квадратурной суммой</i>. Точки <tex>x_i</tex> называются <i>узлами квадратурной формулы</i>. | ||
+ | Разность | ||
+ | |||
+ | :<tex>\Psi _n = \int\limits_a^b f(x)\,dx-\sum_{k=0}^n c_k f(x_k)</tex> | ||
+ | |||
+ | называется <i>погрешностью квадратурной формулы</i>. Погрешность зависит как от расположения узлов, так и от выбора коэффициентов. | ||
== Изложение метода == | == Изложение метода == | ||
- | == | + | === Общие сведения === |
+ | |||
+ | Предположим, что для вычисления интеграла {{eqref|1}} отрезок <tex>[a, b]</tex> разбит на <tex>N</tex> равных отрезков длины | ||
+ | <tex>h = \frac{b-a}N</tex> и на каждом частичном отрезке применяется одна и та же квадратурная формула. Тогда исходный интеграл <tex>I</tex> | ||
+ | заменяется некоторой квадратурной суммой <tex>I_h</tex>, причём возникающая погрешность зависит от шага сетки <tex>h</tex>. | ||
+ | Для некоторых квадратурных формул удаётся получить разложение погрешности <tex>I_h - I</tex> по степеням <tex>h</tex>. Предположим, | ||
+ | что для данной квадратурной суммы <tex>I_h</tex> существует разложение: | ||
+ | |||
+ | {{ eqno | 3}} | ||
+ | :<tex>I_h = I_0 + a_1 h^{\alpha _1} + a_2 h^{\alpha _2} + \ldots + a_m h^{\alpha _m} + O(h^{\alpha _{m+1}})</tex>, | ||
+ | |||
+ | где <tex>0 < \alpha _1 < \alpha _2 < \ldots < \alpha _m < \alpha _{m+1}</tex> и коэффициенты <tex>\{ a_i \} \subset \mathbb{R}</tex> не зависят от <tex>h</tex>. | ||
+ | При этом величины <tex>\{ \alpha _i \} \subset \mathbb{R}</tex> предполагаются известными. | ||
+ | Теперь предположим: | ||
+ | :<tex>I(\frac{h}{r}) = I_0 + a_1\,\frac{h^{\alpha _1}}{r^{\alpha 1}} + a_2\,\frac{h^{\alpha _2}}{r^{\alpha 2}} + \ldots</tex> | ||
+ | |||
+ | Чтобы избавиться от степени <tex>h^{\alpha _1}</tex>, составляющей ошибку (ибо среди всех слагаемых, составляющих ошибку, слагаемое при <tex>h^{\alpha _1}</tex> является наибольшим), вычислим величину <tex>r^{\alpha _1}\,I(\frac{h}r) - I(h)</tex>. Имеем: | ||
+ | |||
+ | :<tex>r^{\alpha _1}\,I(\frac{h}{r}) - I(h) = r^{\alpha _1}\,I_0 - I_0 + a_1 h^{\alpha _1} - a_1 h^{\alpha _1} + a_2\,\frac{h^{\alpha _2}}{r^{\alpha _2 - \alpha _1}} - a_2 h^{\alpha _2} + \ldots</tex> | ||
+ | |||
+ | Отсюда | ||
+ | |||
+ | :<tex>I_1(h) = \frac{r^{\alpha _1}I(\frac{h}r) - I(h)}{r^{\alpha _1} - 1} = I_0 + a_2\,\frac{r^{\alpha _1 - \alpha _2} - 1}{r^{\alpha _1} - 1}\,h^{\alpha _2} + \ldots </tex> | ||
+ | |||
+ | то есть имеем более точное приближение к интегралу <tex>I</tex>. | ||
+ | |||
+ | Таким образом, рекуррентную формулу можно записать в виде: | ||
+ | |||
+ | :<tex>I_{i+1}(h) = I_i(\frac{h}r) + \frac{I_i(\frac{h}r) - I_i(h)}{r^{\alpha _1} - 1} </tex> | ||
+ | |||
+ | Заметим, что <tex>r</tex> - величина, на которую мы делим размер шага при каждом новом вычислении <tex>I</tex>. Разумно положить <tex>r = 2</tex>, т.к. большие значения <tex>r</tex> могут вызвать резкое увеличение количества вычислений. | ||
+ | |||
+ | Для наглядности представим процесс экстраполирования следующей таблицей: | ||
+ | |||
+ | <br> | ||
+ | <tex>\line(1, 0){200}</tex> | ||
+ | <table> | ||
+ | <tr> | ||
+ | <td><tex>I_0^{0}(h)</tex></td> | ||
+ | </tr> | ||
+ | <tr> | ||
+ | <td><tex>I_0^{1}(\frac{h}r)</tex></td> | ||
+ | <td><tex>I_1^{0}(h)</tex></td> | ||
+ | </tr> | ||
+ | <tr> | ||
+ | <td><tex>I_0^{2}(\frac{h}{r^2})</tex></td> | ||
+ | <td><tex>I_1^{1}(\frac{h}r)</tex></td> | ||
+ | <td><tex>I_2^{0}(h)</tex></td> | ||
+ | </tr> | ||
+ | <tr> | ||
+ | <td width = 25%><tex>\vdots</tex></td> | ||
+ | <td width = 25%><tex>\vdots</tex></td> | ||
+ | <td width = 25%><tex>\vdots</tex></td> | ||
+ | <td width = 25%><tex>\ddots</tex></td> | ||
+ | </tr> | ||
+ | </table> | ||
+ | <tex>\line(1,0){200}</tex> | ||
+ | |||
+ | === Метод Рунге === | ||
+ | |||
+ | Пусть существует асимптотическое разложение вида: | ||
+ | |||
+ | <tex>D_N(f) = I_N(f) - I_0 = \alpha _2 h^2 + \alpha _4 h^4 + \alpha _6 h^6 + \ldots</tex>, | ||
+ | |||
+ | где <tex>f(x)</tex> - достаточно гладкая функция, а <tex>hN = b - a</tex>. | ||
+ | |||
+ | Проведем расчеты на двух равномерных сетках с шагами <tex>h_1</tex> и <tex>h_2</tex> соответственно и найдем выражения | ||
+ | <tex>I^{h_1}[f] = I_{N_1}[f]</tex> и <tex>I^{h_2}[f] = I_{N_2}[f]</tex>, <tex>h_1 N_1 = h_2 N_2 = b - a</tex>. Потребуем, | ||
+ | чтобы погрешность для их линейной комбинации: | ||
+ | |||
+ | :<tex>{\tilde D}^{h}(f) = \sigma D^{h_1}(f) + (1 - \sigma)D^{h_2}(f)</tex> | ||
+ | |||
+ | была величиной более высокого порядка по сравнению с <tex>D^{h_1}</tex> и <tex>D^{h_2}</tex>. Если для <tex>D^h = D_N</tex> имеет место формула вида | ||
+ | |||
+ | :<tex> D^{h} = I^{h}[f] - I[f] = \alpha _p h^p + \alpha _q h^q + \ldots,</tex> <tex>q > p</tex>, | ||
+ | |||
+ | то для | ||
+ | |||
+ | :<tex>{\tilde D}^{h} = (\sigma I^{h_1}[f] + (1 - \sigma)I^{h_2}[f]) - I[f]</tex> | ||
+ | |||
+ | получим | ||
+ | |||
+ | :<tex>{\tilde D}^{h}(f) = \alpha _p(\sigma h_1^{p} + (1 - \sigma)h_2^{p}) + \alpha _q (\sigma h_1^{q} + (1 - \sigma)h_2^{q}) + \ldots</tex> | ||
+ | |||
+ | Выберем параметр \sigma из условия <tex>\sigma h_1^{p} + (1 - \sigma)h_2^{p}</tex>: | ||
+ | |||
+ | :<tex>\sigma = \frac{h_2^p}{h_2^p - h_1^p}</tex>. | ||
+ | |||
+ | Имеем | ||
+ | |||
+ | :<tex>{\tilde D}^h(f) = \alpha _q(\sigma h_1^q + (1 - \sigma)h_2^q) + \ldots = O(h^q)</tex>, <tex>h = max(h_1, h_2)</tex>, | ||
+ | |||
+ | причем <tex>\sigma h_1^q + (1 - \sigma)h^q_2 < 0</tex>. Так, если <tex>p = 2, q = 4</tex>, то | ||
+ | <tex>{\tilde D}^{h}(f) = - \alpha _4 h_1^2 h_2^2 + \ldots = O(h^4)</tex>. Таким образом, проведя вычисления на двух сетках | ||
+ | с шагами <tex>h_1</tex> и <tex>h_2 \ne h_1</tex>, мы повысили порядок точности на <tex>2</tex> (на <tex>(q - p)</tex>) для | ||
+ | <tex>{\tilde I} = \sigma I^{h_1} + (1 - \sigma)I^{h_2}</tex>. | ||
+ | |||
+ | === Процесс Эйткена === | ||
+ | |||
+ | Метод расчета на нескольких сетках применяется для повышения порядка точности даже в том случае, когда неизвестен порядок главного члена погрешности. Предположим, что для погрешности имеет место представление | ||
+ | |||
+ | :<tex>D^{h}(f) = \alpha _p h^{p} + \alpha _q h^{q} + \ldots</tex>, <tex>q > p</tex>, | ||
+ | |||
+ | так что | ||
+ | |||
+ | :<tex>I^{h}[f] = I[f] + \alpha _p h^p + \alpha _q h^{q} + \ldots</tex> | ||
+ | |||
+ | Проведем вычисления на трех сетках: <tex>h_1 = h</tex>, <tex>h_2 = \rho h</tex>, <tex>h_3 = \rho^{2} h</tex> (<tex>0 < \rho < 1</tex>). Определим <tex>p</tex>. При этом пренебрегаем членом <tex>O(h^q)</tex>. Образуем отношение | ||
+ | |||
+ | :<tex>A = \frac{I^{h_1}[f] - I^{h_2}[f]}{I^{h_2}[f] - I^{h_3}[f]} \approx \frac{h^p_1 - h^p_2}{h^p_2 - h^p_3} = \frac{1 - \rho^p}{\rho^p(1-\rho^p)} = (\frac{1}{\rho})^p</tex> | ||
+ | |||
+ | и найдем | ||
+ | |||
+ | :<tex>p \approx \frac{\ln A}{\ln(\frac{1}{\rho})}</tex>. | ||
+ | |||
+ | Зная приближенное значение <tex>p</tex>, можно методом Рунге повысить порядок точности. Для этого образуем комбинацию | ||
+ | <tex>{\tilde I}^h = \sigma I^{h_1} + (1 - \sigma)I^{h_2}</tex> и выберем <tex>\sigma</tex> так, чтобы <tex> \sigma h_1^p+(1 - \sigma)h_2^p = (\sigma + (1 - \sigma)\rho^p)h^p = 0 </tex>, то есть <tex>\sigma = \frac{\rho^p}{\rho^p - 1} = \frac{1}{1 - A}</tex>. | ||
+ | Тогда для погрешности <tex>{\tilde D}^h = {\tilde I}^h - I</tex> получаем | ||
+ | |||
+ | :<tex> {\tilde D}^h(f) = O(h^q) </tex>. | ||
+ | |||
+ | |||
== Числовой пример == | == Числовой пример == | ||
+ | |||
+ | Найдем с помощью квадратурной формулы трапеций приближенное значение интеграла, применив экстраполяцию Ричардсона (данный метод называется <i>методом Ромберга</i>): | ||
+ | |||
+ | :<tex>\int^{5}_{1} \ln x\, dx = x \ln x - x |_1^{5} = 5\ln 5 - 4</tex> | ||
+ | |||
+ | [[Изображение:Extrapol.png|thumb|300px]] | ||
+ | В нижеследующей таблице представлены результаты работы [[Media:Richardson.zip|программы]]: | ||
+ | {|class="standard" | ||
+ | !r | ||
+ | !Исходная формула | ||
+ | !Экстраполированная формула | ||
+ | !Точное значение | ||
+ | !Погрешность вычислений | ||
+ | !Погрешность формулы | ||
+ | |- | ||
+ | !2 | ||
+ | |3.98277278 | ||
+ | |4.04665506 | ||
+ | |4.04718956 | ||
+ | |0.0005345 | ||
+ | |0.00275556 | ||
+ | |- | ||
+ | !4 | ||
+ | |4.03068449 | ||
+ | |4.04714980 | ||
+ | |4.04718956 | ||
+ | |0.00003976 | ||
+ | |0.00017222 | ||
+ | |- | ||
+ | !8 | ||
+ | |4.04303347 | ||
+ | |4.04718692 | ||
+ | |4.04718956 | ||
+ | |0.00000264 | ||
+ | |0.00001076 | ||
+ | |- | ||
+ | !16 | ||
+ | |4.04614856 | ||
+ | |4.04718939 | ||
+ | |4.04718956 | ||
+ | |0.00000017 | ||
+ | |0.00000067 | ||
+ | |- | ||
+ | !32 | ||
+ | |4.04692918 | ||
+ | |4.04718955 | ||
+ | |4.04718956 | ||
+ | |0.00000001 | ||
+ | |0.00000004 | ||
+ | |- | ||
+ | !64 | ||
+ | |4.04712446 | ||
+ | |4.04718956 | ||
+ | |4.04718956 | ||
+ | |0 | ||
+ | |0 | ||
+ | |- | ||
+ | !20384 | ||
+ | |4.04718956 | ||
+ | | | ||
+ | | | ||
+ | | | ||
+ | | | ||
+ | |} | ||
+ | |||
+ | Здесь <tex>r</tex> - коэффициент измельчения шага <tex>h</tex>. Исходная величина шага <tex>h = 2</tex>. | ||
+ | |||
+ | На иллюстрации черная сплошная линия - исходная формула, красная пунктирная - экстраполированная. | ||
+ | |||
+ | Как мы видим, разница между экстраполированными и неэкстраполированными результатами значительна. Уже при величине шага в <tex>h = \frac{1}{64}</tex> мы можем найти значение интеграла с точностью <tex>10^{-8}</tex>, тогда как в исходной формуле нам для достижения такой точности пришлось бы задать величину шага <tex>h = \frac{1}{10192}</tex>. | ||
== Рекомендации программисту == | == Рекомендации программисту == | ||
+ | |||
+ | Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии. | ||
+ | Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца. | ||
== Заключение == | == Заключение == | ||
- | + | Мы получили более точные результаты при меньшем количестве вычислений функции, чем в исходном методе. Избавившись от степени, составляющей ошибку, мы пришли к результату, который в противном случае был бы недостижим без значительного уменьшения размера шага. Таким образом, мы преобразовали весьма посредственный алгоритм вычисления интеграла в довольно эффективный. | |
- | + | == Список литературы == | |
+ | * ''А.А.Самарский, А.В.Гулин.'' Численные методы М.: Наука, 1989. | ||
+ | * ''А.А.Самарский.'' Введение в численные методы М.: Наука, 1982. | ||
+ | * [http://ocw.mit.edu/NR/rdonlyres/Mathematics/18-304Spring-2006/D69D4111-2100-443D-A75D-5D5BE7B4FAA2/0/xtrpltn_liu_xpnd.pdf Fundamental Methods of Numerical Extrapolation With Applications], mit.edu | ||
+ | == См. также == | ||
+ | * [[Методы прямоугольников и трапеций]] | ||
+ | * [[Методы парабол (Симпсона) и более высоких степеней (Ньютона - Котеса)]] | ||
[[Категория:Численное интегрирование]] | [[Категория:Численное интегрирование]] | ||
+ | [[Категория:Учебные задачи]] |
Текущая версия
Содержание |
Введение
Постановка математической задачи
Задача численного интегрирования состоит в приближенном нахождении значения интеграла
где
- заданная и интегрируемая на
функция. В качестве приближенного значения рассматривается число
где - числовые коэффициенты и
- точки отрезка
,
.
Приближенное равенство
называется квадратурной формулой, а сумма вида (2) - квадратурной суммой. Точки называются узлами квадратурной формулы.
Разность
называется погрешностью квадратурной формулы. Погрешность зависит как от расположения узлов, так и от выбора коэффициентов.
Изложение метода
Общие сведения
Предположим, что для вычисления интеграла (1) отрезок разбит на
равных отрезков длины
и на каждом частичном отрезке применяется одна и та же квадратурная формула. Тогда исходный интеграл
заменяется некоторой квадратурной суммой
, причём возникающая погрешность зависит от шага сетки
.
Для некоторых квадратурных формул удаётся получить разложение погрешности
по степеням
. Предположим,
что для данной квадратурной суммы
существует разложение:
,
где и коэффициенты
не зависят от
.
При этом величины
предполагаются известными.
Теперь предположим:
Чтобы избавиться от степени , составляющей ошибку (ибо среди всех слагаемых, составляющих ошибку, слагаемое при
является наибольшим), вычислим величину
. Имеем:
Отсюда
то есть имеем более точное приближение к интегралу .
Таким образом, рекуррентную формулу можно записать в виде:
Заметим, что - величина, на которую мы делим размер шага при каждом новом вычислении
. Разумно положить
, т.к. большие значения
могут вызвать резкое увеличение количества вычислений.
Для наглядности представим процесс экстраполирования следующей таблицей:
Метод Рунге
Пусть существует асимптотическое разложение вида:
,
где - достаточно гладкая функция, а
.
Проведем расчеты на двух равномерных сетках с шагами и
соответственно и найдем выражения
и
,
. Потребуем,
чтобы погрешность для их линейной комбинации:
была величиной более высокого порядка по сравнению с и
. Если для
имеет место формула вида
,
то для
получим
Выберем параметр \sigma из условия :
.
Имеем
,
,
причем . Так, если
, то
. Таким образом, проведя вычисления на двух сетках
с шагами
и
, мы повысили порядок точности на
(на
) для
.
Процесс Эйткена
Метод расчета на нескольких сетках применяется для повышения порядка точности даже в том случае, когда неизвестен порядок главного члена погрешности. Предположим, что для погрешности имеет место представление
,
,
так что
Проведем вычисления на трех сетках: ,
,
(
). Определим
. При этом пренебрегаем членом
. Образуем отношение
и найдем
.
Зная приближенное значение , можно методом Рунге повысить порядок точности. Для этого образуем комбинацию
и выберем
так, чтобы
, то есть
.
Тогда для погрешности
получаем
.
Числовой пример
Найдем с помощью квадратурной формулы трапеций приближенное значение интеграла, применив экстраполяцию Ричардсона (данный метод называется методом Ромберга):
В нижеследующей таблице представлены результаты работы программы:
r | Исходная формула | Экстраполированная формула | Точное значение | Погрешность вычислений | Погрешность формулы |
---|---|---|---|---|---|
2 | 3.98277278 | 4.04665506 | 4.04718956 | 0.0005345 | 0.00275556 |
4 | 4.03068449 | 4.04714980 | 4.04718956 | 0.00003976 | 0.00017222 |
8 | 4.04303347 | 4.04718692 | 4.04718956 | 0.00000264 | 0.00001076 |
16 | 4.04614856 | 4.04718939 | 4.04718956 | 0.00000017 | 0.00000067 |
32 | 4.04692918 | 4.04718955 | 4.04718956 | 0.00000001 | 0.00000004 |
64 | 4.04712446 | 4.04718956 | 4.04718956 | 0 | 0 |
20384 | 4.04718956 |
Здесь - коэффициент измельчения шага
. Исходная величина шага
.
На иллюстрации черная сплошная линия - исходная формула, красная пунктирная - экстраполированная.
Как мы видим, разница между экстраполированными и неэкстраполированными результатами значительна. Уже при величине шага в мы можем найти значение интеграла с точностью
, тогда как в исходной формуле нам для достижения такой точности пришлось бы задать величину шага
.
Рекомендации программисту
Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии. Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца.
Заключение
Мы получили более точные результаты при меньшем количестве вычислений функции, чем в исходном методе. Избавившись от степени, составляющей ошибку, мы пришли к результату, который в противном случае был бы недостижим без значительного уменьшения размера шага. Таким образом, мы преобразовали весьма посредственный алгоритм вычисления интеграла в довольно эффективный.
Список литературы
- А.А.Самарский, А.В.Гулин. Численные методы М.: Наука, 1989.
- А.А.Самарский. Введение в численные методы М.: Наука, 1982.
- Fundamental Methods of Numerical Extrapolation With Applications, mit.edu