Участник:Gukov/Песочница

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

(Различия между версиями)
Перейти к: навигация, поиск
(Постановка математической задачи)
(Случай почти равных корней)
 
(29 промежуточных версий не показаны.)
Строка 1: Строка 1:
== Введение ==
== Введение ==
-
=== Постановка математической задачи ===
+
=== Метод Ньютона-Рафсона ===
 +
{{main|Метод касательных (Ньютона-Рафсона)}}
-
Задача численного интегрирования состоит в приближенном нахождении значения интеграла
+
Пусть имеется некоторая функция <tex>F(x)</tex> и необходимо найти такие значения аргумента <tex>x</tex>, для которых
-
{{ eqno | 1}}
+
{{eqno|1}}
-
:<tex>I = \int\limits_a^b f(x)\,dx,</tex>
+
:<tex>F(x)=0</tex>
-
где
+
Перепишем {{eqref|1}} в виде
-
<tex>
+
{{eqno|1.1}}
-
f(x)
+
:<tex>x = f(x)</tex>
-
</tex> - заданная и интегрируемая на <tex> [a, b] </tex> функция. В качестве приближенного значения рассматривается число
+
-
{{ eqno | 2}}
+
и запишем <tex>n+1</tex>-ое приближения корня {{eqref|1.1}}, при этом делая поправку <tex>\alpha</tex> к очередному значению <tex>x_n</tex>
-
:<tex>I_n=\sum_{i=0}^n c_k f(x_k),</tex>
+
{{eqno|2}}
 +
:<tex>x_{n+1} = x_n + \alpha \Delta x_n,</tex>
 +
где <tex>\Delta x_n = f(x_n) - x_n,</tex> и положим <tex>\alpha = \frac{1}{1-f'(x_n)}</tex>.
-
где <tex>c_k</tex> - числовые коэффициенты и <tex>x_k</tex> - точки отрезка <tex>[a,b]</tex>, <tex> k = 0, 1, \ldots, n </tex>.
+
Тогда {{eqref|2}} перепишется в виде
-
Приближенное равенство
+
{{eqno|3}}
 +
:<tex>x_{n+1} = \frac{f(x_n) - x_n f'(x_n)}{1 - f'(x_n)}</tex>
-
:<tex>\int\limits_a^b f(x)\,dx=\sum_{k=0}^n c_k f(x_k)</tex>
+
Нетрудно видеть, что {{eqref|3}} эквивалентно простому [[Метод последовательных приближений|методу последовательных приближений]]
 +
:<tex>x_{n+1} = g(x_n),</tex>
 +
где <tex>g(x) = \frac{f(x) - x f'(x)}{1 - f'(x)}</tex>.
-
называется <i>квадратурной формулой</i>, а сумма вида {{eqref|2}} - <i>квадартурной суммой</i>. Точки <tex>x_i</tex> называются <i>узлами квадратурной формулы</i>.
+
Вспомним также, что если <tex>|g'(x)| < 1</tex>, то метод сходится. Имеем
-
Разность
+
:<tex>g'(x) = \frac{f''(x)[f(x) - x]}{[1 - f'(x)]^2}.</tex>
-
:<tex>\Psi _n = \int\limits_a^b f(x)\,dx-\sum_{k=0}^n c_k f(x_k)</tex>
+
Но так как справедливо соотношение {{eqref|1.1}}, то для <tex>x,</tex> достаточно близких к решению {{eqref|1.1}}, выражение в скобках в числителе дроби становится малым. Поэтому итерационный метод, описываемый формулой {{eqref|3}} сходится, если
 +
* 1. Начальное приближение <tex>x_0</tex> выбрано достаточно близко к решению <tex>x = f(x)</tex>.
 +
* 2. Производная <tex>f''(x)</tex> не становится слишком большой.
 +
* 3. Производная <tex>f'(x)</tex> не слишком близка к 1.
-
называется <i>погрешностью квадратурной формулы</i>. Погрешность зависит как от расположения узлов, так и от выбора коэффициентов.
+
Это и есть метод Ньютона-Рафсона. Обычно его записывают в виде
 +
{{eqno|4}}
 +
:<tex>x_{n+1} = x_{n} - \frac{F(x_n)}{F'(x_n)}</tex>,
 +
 
 +
где <tex>F(x) = f(x) - x = 0</tex>
 +
 
 +
Таким образом, мы вернулись к уравнению в форме {{eqref|1}}, и условия сходимости принимают следующий вид
 +
* 1. Начальное приближение <tex>x_0</tex> выбрано достаточно близко к корню уравнения <tex>F(x) = 0</tex>.
 +
* 2. Производная <tex>F''(x)</tex> не становится очень большой.
 +
* 3. Производная <tex>F'(x)</tex> не слишком близка к 0.
 +
 
 +
=== Случай почти равных корней ===
 +
[[Изображение:macon.png|thumb|200px|Рисунок 1]]
 +
Условие 3 сходимости метода Ньютона-Рафсона означает, что никакие два корня не находятся слишком близко один к другому. Соответствующая ситуация представлена на рисунке 1 (масштаб сильно увеличен). Заметим, что производная <tex>f'(x)</tex> близка
 +
к 1 при x, равном обоим значениям корней, <tex>a_1</tex> и <tex>a_2</tex>. Более того, на основании теоремы Лагранжа, можно утверждать, что <tex>f'(x) = 1</tex> где-то между <tex>a_1</tex> и <tex>a_2</tex>.
 +
 
 +
Рассмотрим, что случится, если принять <tex>x_0</tex> в качестве исходного значения для корня <tex>a_1</tex>. Касательная,
 +
проведенная через точку <tex>C</tex>, пересечет прямую <tex>y=x</tex> в точке <tex>A</tex>, и следующее приближение будет равно
 +
<tex>x_1</tex>. Касательная, проведенная через точку <tex>B</tex>, пересекает прямую в точке <tex>D</tex>, и в качестве следующего приближеня получается снова <tex>x_0</tex>. Итерационный процесс, таким образом, осциллирует между <tex>x_0</tex> и <tex>x_1</tex>
 +
до бесконечности, не сходясь ни к одному значению корня. Иначе говоря, не удается <i>отделить</i> эти два корня, потому что они расположены слишком близко.
 +
 
 +
Поэтому необходимо начальное приближение, достаточно близкое к искомому значению корня. Трудности возникают потому, что
 +
вычисление знаменателя в формуле {{eqref|3}} включает в себя вычитание двух почти равных чисел, что приводит к понижению точности.
== Изложение метода ==
== Изложение метода ==
 +
 +
=== Поиск начального приближения ===
 +
 +
[[Изображение:macon2.PNG|thumb|200px|Рисунок 2]]
 +
Сначала находится значение x, при котором <tex>f'(x)=1</tex>, то есть решается уравнение
 +
:<tex>x=x+f'(x)-1</tex>
 +
Пусть решением этого уравнения будет некоторое <tex>x = x_0</tex>. Эта точка расположена между двумя корнями, <tex>a_1</tex> и <tex>a_2</tex>. Чтобы получить начальное приближение для решения уравнения, предположим, что <tex>x_0</tex> лежит посредине между
 +
<tex>a_1</tex> и <tex>a_2</tex> (рисунок 2). Другими словами, мы предполагаем, что <tex>x_0 + d</tex> и <tex>x_0 - d</tex> являются корнями уравнения {{eqref|1.1}}. Разлагая <tex>f(x)</tex> в ряд Тейлора в окрестности точки <tex>x_0</tex> и принимая во внимание, что <tex>f'(x)=1</tex>, получаем
 +
:<tex>f(x) = f(x_0) + (x - x_0) + \frac{1}{2} f''(x_0) (x - x_0)^2 + \ldots</tex>
 +
 +
Ограничим ряд тремя членами. Подставляя <tex>x + d</tex> вместо <tex>x</tex>, имеем
 +
:<tex>f(x_0+d) = f(x_0) + d + \frac{1}{2} f''(x_0)(d^2)</tex>
 +
 +
Но по условию
 +
:<tex>f(x_0+d) = x_0+d,</tex>
 +
 +
поэтому, решая эти уравнения относительно <tex>d</tex>, получаем
 +
:<tex>d=\sqrt{\frac{2(x_0 - f(x_0))}{f''(x_0)}}</tex>
 +
 +
Таким образом, процесс решения сводится к следующему. Если дано уравнение с почти равными корнями, то, определив приблизительное
 +
местонахождение этих корней, необходиом решить уравнение
 +
:<tex>x = x + f'(x) - 1</tex>
 +
 +
и определить значение <tex>x_0</tex>. Для решения этого уравнения можно применить, например, метод Ньютона-Рафсона.
 +
Найдя значение <tex>x_0</tex>, можно определить значение <tex>d</tex>. И, наконец, значения <tex>x_0-d</tex> и <tex>x_0+d</tex>
 +
используются в качестве начальных приближений для определения соответственно <tex>a_1</tex> и <tex>a_2</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>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<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>, &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <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>,&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<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]]
 +
В нижеследующей таблице представлены результаты работы программы:
 +
{|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>.
== Рекомендации программисту ==
== Рекомендации программисту ==
 +
 +
Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии.
 +
Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца.
== Заключение ==
== Заключение ==
-
== Список литературы ==
+
Мы получили более точные результаты при меньшем количестве вычислений функции, чем в исходном методе. Избавившись от степени, составляющей ошибку, мы пришли к результату, который в противном случае был бы недостижим без значительного уменьшения размера шага. Таким образом, мы преобразовали весьма посредственный алгоритм вычисления интегралов в довольно эффективный.
 +
== Список литературы ==
 +
* ''А.А.Самарский, А.В.Гулин.''&nbsp; Численные методы М.: Наука, 1989.
 +
* [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
{{Stub}}
{{Stub}}

Текущая версия

Содержание

Введение

Метод Ньютона-Рафсона

Пусть имеется некоторая функция F(x) и необходимо найти такие значения аргумента x, для которых

(1)
F(x)=0

Перепишем (1) в виде

(1.1)
x = f(x)

и запишем n+1-ое приближения корня (1.1), при этом делая поправку \alpha к очередному значению x_n

(2)
x_{n+1} = x_n + \alpha \Delta x_n,

где \Delta x_n = f(x_n) - x_n, и положим \alpha = \frac{1}{1-f'(x_n)}.

Тогда (2) перепишется в виде

(3)
x_{n+1} = \frac{f(x_n) - x_n f'(x_n)}{1 - f'(x_n)}

Нетрудно видеть, что (3) эквивалентно простому методу последовательных приближений

x_{n+1} = g(x_n),

где g(x) = \frac{f(x) - x f'(x)}{1 - f'(x)}.

Вспомним также, что если |g'(x)| < 1, то метод сходится. Имеем

g'(x) = \frac{f''(x)[f(x) - x]}{[1 - f'(x)]^2}.

Но так как справедливо соотношение (1.1), то для x, достаточно близких к решению (1.1), выражение в скобках в числителе дроби становится малым. Поэтому итерационный метод, описываемый формулой (3) сходится, если

  • 1. Начальное приближение x_0 выбрано достаточно близко к решению x = f(x).
  • 2. Производная f''(x) не становится слишком большой.
  • 3. Производная f'(x) не слишком близка к 1.

Это и есть метод Ньютона-Рафсона. Обычно его записывают в виде

(4)
x_{n+1} = x_{n} - \frac{F(x_n)}{F'(x_n)},

где F(x) = f(x) - x = 0

Таким образом, мы вернулись к уравнению в форме (1), и условия сходимости принимают следующий вид

  • 1. Начальное приближение x_0 выбрано достаточно близко к корню уравнения F(x) = 0.
  • 2. Производная F''(x) не становится очень большой.
  • 3. Производная F'(x) не слишком близка к 0.

Случай почти равных корней

Рисунок 1
Рисунок 1

Условие 3 сходимости метода Ньютона-Рафсона означает, что никакие два корня не находятся слишком близко один к другому. Соответствующая ситуация представлена на рисунке 1 (масштаб сильно увеличен). Заметим, что производная f'(x) близка к 1 при x, равном обоим значениям корней, a_1 и a_2. Более того, на основании теоремы Лагранжа, можно утверждать, что f'(x) = 1 где-то между a_1 и a_2.

Рассмотрим, что случится, если принять x_0 в качестве исходного значения для корня a_1. Касательная, проведенная через точку C, пересечет прямую y=x в точке A, и следующее приближение будет равно x_1. Касательная, проведенная через точку B, пересекает прямую в точке D, и в качестве следующего приближеня получается снова x_0. Итерационный процесс, таким образом, осциллирует между x_0 и x_1 до бесконечности, не сходясь ни к одному значению корня. Иначе говоря, не удается отделить эти два корня, потому что они расположены слишком близко.

Поэтому необходимо начальное приближение, достаточно близкое к искомому значению корня. Трудности возникают потому, что вычисление знаменателя в формуле (3) включает в себя вычитание двух почти равных чисел, что приводит к понижению точности.

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

Поиск начального приближения

Рисунок 2
Рисунок 2

Сначала находится значение x, при котором f'(x)=1, то есть решается уравнение

x=x+f'(x)-1

Пусть решением этого уравнения будет некоторое x = x_0. Эта точка расположена между двумя корнями, a_1 и a_2. Чтобы получить начальное приближение для решения уравнения, предположим, что x_0 лежит посредине между a_1 и a_2 (рисунок 2). Другими словами, мы предполагаем, что x_0 + d и x_0 - d являются корнями уравнения (1.1). Разлагая f(x) в ряд Тейлора в окрестности точки x_0 и принимая во внимание, что f'(x)=1, получаем

f(x) = f(x_0) + (x - x_0) + \frac{1}{2} f''(x_0) (x - x_0)^2 + \ldots

Ограничим ряд тремя членами. Подставляя x + d вместо x, имеем

f(x_0+d) = f(x_0) + d + \frac{1}{2} f''(x_0)(d^2)

Но по условию

f(x_0+d) = x_0+d,

поэтому, решая эти уравнения относительно d, получаем

d=\sqrt{\frac{2(x_0 - f(x_0))}{f''(x_0)}}

Таким образом, процесс решения сводится к следующему. Если дано уравнение с почти равными корнями, то, определив приблизительное местонахождение этих корней, необходиом решить уравнение

x = x + f'(x) - 1

и определить значение x_0. Для решения этого уравнения можно применить, например, метод Ньютона-Рафсона. Найдя значение x_0, можно определить значение d. И, наконец, значения x_0-d и x_0+d используются в качестве начальных приближений для определения соответственно a_1 и a_2.

Метод Рунге

Пусть существует асимптотическое разложение вида:

D_N(f) = I_N(f) - I_0 = \alpha _2 h^2 + \alpha _4 h^4 + \alpha _6 h^6 + \ldots,

где f(x) - достаточно гладкая функция, а hN = b - a.

Проведем расчеты на двух равномерных сетках с шагами h_1 и h_2 соответственно и найдем выражения I^{h_1}[f] = I_{N_1}[f] и I^{h_2}[f] = I_{N_2}[f], h_1 N_1 = h_2 N_2 = b - a. Потребуем, чтобы погрешность для их линейной комбинации:

{\tilde D}^{h}(f) = \sigma D^{h_1}(f) + (1 - \sigma)D^{h_2}(f)

была величиной более высокого порядка по сравнению с D^{h_1} и D^{h_2}. Если для D^h = D_N имеет место формула вида

 D^{h} = I^{h}[f] - I[f] = \alpha _p h^p + \alpha _q h^q + \ldots,     q > p,

то для

{\tilde D}^{h} = (\sigma I^{h_1}[f] + (1 - \sigma)I^{h_2}[f]) - I[f]

получим

{\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

Выберем параметр \sigma из условия \sigma h_1^{p} + (1 - \sigma)h_2^{p}:

\sigma = \frac{h_2^p}{h_2^p - h_1^p}.

Имеем

{\tilde D}^h(f) = \alpha _q(\sigma h_1^q + (1 - \sigma)h_2^q) + \ldots = O(h^q),       h = max(h_1, h_2),

причем \sigma h_1^q + (1 - \sigma)h^q_2 < 0. Так, если p = 2, q = 4, то {\tilde D}^{h}(f) = - \alpha _4 h_1^2 h_2^2 + \ldots = O(h^4). Таким образом, проведя вычисления на двух сетках с шагами h_1 и h_2 \ne h_1, мы повысили порядок точности на 2 (на (q - p)) для {\tilde I} = \sigma I^{h_1} + (1 - \sigma)I^{h_2}.

Процесс Эйткена

Метод расчета на нескольких сетках применяется для повышения порядка точности даже в том случае, когда неизвестен порядок главного члена погрешности. Предположим, чт для погрешности имеет место представление

D^{h}(f) = \alpha _p h^{p} + \alpha _q h^{q} + \ldots,     q > p,

так что

I^{h}[f] = I[f] + \alpha _p h^p + \alpha _q h^{q} + \ldots

Проведем вычисления на трех сетках: h_1 = h, h_2 = \rho h, h_3 = \rho^{2} h (0 < \rho < 1). Определим p. При этом пренебрегаем членом O(h^q). Образуем отношение

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

и найдем

p \approx \frac{\ln A}{\ln(\frac{1}{\rho})}.

Зная приближенное значение p, можно методом методом Рунге повысить порядок точности. Для этого образуем комбинацию {\tilde I}^h = \sigma I^{h_1} + (1 - \sigma)I^{h_2} и выберем \sigma так, чтобы  \sigma h_1^p+(1 - \sigma)h_2^p = (\sigma + (1 - \sigma)\rho^p)h^p = 0 , то есть \sigma = \frac{\rho^p}{\rho^p - 1} = \frac{1}{1 - A}. Тогда для погрешности {\tilde D}^h = {\tilde I}^h - I получаем

 {\tilde D}^h(f) = O(h^q) .

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

Найдем с помощью квадратурной формулы трапеций приближенное значение интеграла, применив экстраполяцию Ричардсона (данный метод называется методом Ромберга):

\int^{5}_{1} \ln x\, dx = x \ln x - x |_1^{5} = 5\ln 5 - 4

В нижеследующей таблице представлены результаты работы программы:

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

Здесь r - коэффициент измельчения шага h. Исходная величина шага h = 2.

На иллюстрации черная сплошная линия - исходная формула, красная пунктирная - экстраполированная.

Как мы видим, разница между экстраполированными и неэкстраполированными результатами значительна. Уже при величине шага в h = \frac{1}{64} мы можем найти значение интеграла с точностью 10^{-8}, тогда как в исходной формуле нам для достижения такой точности пришлось бы задать величину шага h = \frac{1}{10192}.

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

Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии. Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца.

Заключение

Мы получили более точные результаты при меньшем количестве вычислений функции, чем в исходном методе. Избавившись от степени, составляющей ошибку, мы пришли к результату, который в противном случае был бы недостижим без значительного уменьшения размера шага. Таким образом, мы преобразовали весьма посредственный алгоритм вычисления интегралов в довольно эффективный.

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

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