Экстраполяция Ричардсона - Richardson extrapolation

В численный анализ, Экстраполяция Ричардсона это ускорение последовательности метод, используемый для улучшения скорость сходимости из последовательность оценок некоторой стоимости ${ displaystyle A ^ { ast} = lim _ {h to 0} A (h)}$ . По сути, с учетом стоимости ${ displaystyle A (h)}$ для нескольких значений ${ displaystyle h}$ , мы можем оценить ${ displaystyle A ^ { ast}}$ путем экстраполяции оценок на ${ displaystyle h = 0}$ . Он назван в честь Льюис Фрай Ричардсон, который представил технику в начале 20 века.^[1]^[2] По словам Биркофф и Рота, «его полезность для практических вычислений трудно переоценить».^[3]

Практические применения экстраполяции Ричардсона включают: Интеграция Ромберга, который применяет экстраполяцию Ричардсона к правило трапеции, а Алгоритм Булирша – Стоера для решения обыкновенных дифференциальных уравнений.

Пример экстраполяции Ричардсона

Предположим, что мы хотим приблизить ${ displaystyle A ^ {*}}$ , и у нас есть метод ${ displaystyle A (h)}$ это зависит от небольшого параметра ${ displaystyle h}$ таким образом, что

${ displaystyle A (h) = A ^ { ast} + Ch ^ {n} + O (h ^ {n + 1}).}$

Определим новую функцию

${ Displaystyle R (ч, т): = { гидроразрыва {т ^ {п} А (ч / т) -А (ч)} {т ^ {п} -1}}}$ где ${ displaystyle h}$ и ${ displaystyle { frac {h} {t}}}$ - это два разных размера шага.

потом

${ Displaystyle R (час, t) = { гидроразрыва {t ^ {n} (A ^ {*} + C left ({ frac {h} {t}} right) ^ {n} + O ( h ^ {n + 1})) - (A ^ {*} + Ch ^ {n} + O (h ^ {n + 1}))} {t ^ {n} -1}} = A ^ {* } + O (h ^ {n + 1}).}$

${ Displaystyle R (ч, т)}$ называется Ричардсон экстраполяция из А(час), и имеет оценку ошибки более высокого порядка ${ Displaystyle О (ч ^ {п + 1})}$ в сравнении с ${ displaystyle A (h)}$ .

Очень часто гораздо проще получить заданную точность, используя R (ч) скорее, чем А (ч ') с гораздо меньшим час' . куда А (ч ') может вызвать проблемы из-за ограниченной точности (ошибки округления ) и / или из-за увеличения количество расчетов необходимо (см. примеры ниже).

Общая формула

Позволять ${ displaystyle A (h)}$ быть приближением ${ displaystyle A ^ {*}}$ (точное значение), которое зависит от положительного размера шага час с ошибка формула формы

{ displaystyle A ^ {*} = A (h) + a_ {0} h ^ {k_ {0}} + a_ {1} h ^ {k_ {1}} + a_ {2} h ^ {k_ {2 }} + cdots}

где а_я неизвестные константы и k_я - известные константы такие, что час^k_я > час^{k_{я + 1}}.

k₀ - поведение размера шага ведущего порядка Ошибка усечения так как ${ displaystyle A ^ {*} = A (h) + O (h ^ {k_ {0}})}$

Искомое точное значение можно определить как

{ displaystyle A ^ {*} = A (h) + a_ {0} h ^ {k_ {0}} + a_ {1} h ^ {k_ {1}} + a_ {2} h ^ {k_ {2 }} + cdots}

который можно упростить с помощью Обозначение Big O быть

{ displaystyle A ^ {*} = A (h) + a_ {0} h ^ {k_ {0}} + O (h ^ {k_ {1}}). , !}

Использование размеров шага час и ${ displaystyle h / t}$ для некоторой постоянной т, две формулы для А находятся:

{ displaystyle A ^ {*} = A (h) + a_ {0} h ^ {k_ {0}} + O (h ^ {k_ {1}}) , !}

{ displaystyle A ^ {*} = A ! left ({ frac {h} {t}} right) + a_ {0} left ({ frac {h} {t}} right) ^ {k_ {0}} + O (h ^ {k_ {1}}).}

Умножая второе уравнение на т^k₀ и вычитание первого уравнения дает

{ displaystyle (t ^ {k_ {0}} - 1) A ^ {*} = t ^ {k_ {0}} A left ({ frac {h} {t}} right) -A (h ) + (t ^ {k_ {0}} - 1) O (h ^ {k_ {1}})}

который может быть решен для ${ displaystyle A ^ {*}}$ давать

{ displaystyle A ^ {*} = { frac {t ^ {k_ {0}} A left ({ frac {h} {t}} right) -A (h)} {t ^ {k_ { 0}} - 1}} + O (h ^ {k_ {1}}).}

Следовательно, используя ${ displaystyle R (h, t) = { frac {t ^ {k_ {0}} A left ({ frac {h} {t}} right) -A (h)} {t ^ {k_ {0}} - 1}}}$ ошибка усечения была уменьшена до ${ Displaystyle О (ч ^ {к_ {1}})}$ . Это в отличие от ${ displaystyle A (h)}$ где ошибка усечения является ${ Displaystyle О (ч ^ {к_ {0}})}$ для того же размера шага ${ displaystyle h}$

Благодаря этому процессу мы достигли лучшего приближения к А путем вычитания наибольшего члена ошибки, которая была О(час^k₀). Этот процесс можно повторить, чтобы удалить больше ошибок, чтобы получить еще лучшие приближения.

Генерал отношение повторения начиная с ${ displaystyle A_ {0} = A (h)}$ можно определить для приближений как

{ Displaystyle A_ {я + 1} (h) = { frac {t ^ {k_ {i}} A_ {i} left ({ frac {h} {t}} right) -A_ {i} (h)} {t ^ {k_ {i}} - 1}}}

где ${ displaystyle k_ {я + 1}}$ удовлетворяет

{ displaystyle A ^ {*} = A_ {я + 1} (h) + O (h ^ {k_ {i + 1}})}

.

Экстраполяцию Ричардсона можно рассматривать как линейную преобразование последовательности.

Кроме того, общая формула может использоваться для оценки k₀ (поведение размера шага ведущего порядка Ошибка усечения ), когда ни его значение, ни А^* (точное значение) известно априори. Такой метод может быть полезен для количественной оценки неизвестного скорость сходимости. Учитывая приближения А от трех различных размеров шага час, ч / т, и ч / с, точное отношение

{ displaystyle A ^ {*} = { frac {t ^ {k_ {0}} A left ({ frac {h} {t}} right) -A (h)} {t ^ {k_ { 0}} - 1}} + O (h ^ {k_ {1}}) = { frac {s ^ {k_ {0}} A left ({ frac {h} {s}} right) - A (h)} {s ^ {k_ {0}} - 1}} + O (h ^ {k_ {1}})}

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

{ displaystyle A left ({ frac {h} {t}} right) + { frac {A left ({ frac {h} {t}} right) -A (h)} {t ^ {k_ {0}} - 1}} приблизительно A left ({ frac {h} {s}} right) + { frac {A left ({ frac {h} {s}} справа) -A (h)} {s ^ {k_ {0}} - 1}}}

которую можно решить численно для оценки k₀ для произвольного выбора час, s и т.

Пример кода псевдокода для экстраполяции Ричардсона

Следующий псевдокод в стиле MATLAB демонстрирует экстраполяцию Ричардсона, чтобы помочь решить ODE ${ Displaystyle у '(т) = - у ^ {2}}$ , ${ Displaystyle у (0) = 1}$ с Трапециевидный метод. В этом примере мы уменьшаем размер шага вдвое. ${ displaystyle h}$ каждую итерацию, поэтому в приведенном выше обсуждении ${ displaystyle t = 2}$ . Ошибка трапециевидного метода может быть выражена в терминах нечетных степеней, так что погрешность на нескольких шагах может быть выражена в четных степенях; это заставляет нас поднять ${ displaystyle t}$ ко второй власти и взять власть ${ Displaystyle 4 = 2 ^ {2} = т ^ {2}}$ в псевдокоде. Мы хотим найти значение ${ Displaystyle у (5)}$ , которая имеет точное решение ${ displaystyle { frac {1} {5 + 1}} = { frac {1} {6}} = 0,1666 ...}$ так как точное решение ОДУ есть ${ Displaystyle у (т) = { гидроразрыва {1} {1 + т}}}$ . Этот псевдокод предполагает, что функция с именем Трапециевидный (f, tStart, tEnd, h, y0) существует, который пытается вычислить у (конец) путем выполнения трапециевидного метода над функцией ж, с отправной точкой y0 и tStart и размер шага час.

Обратите внимание, что слишком маленький начальный размер шага может потенциально внести ошибку в окончательное решение. Хотя существуют методы, предназначенные для выбора наилучшего начального размера шага, один из вариантов - начать с большого размера шага, а затем позволить экстраполяции Ричардсона уменьшать размер шага на каждой итерации, пока ошибка не достигнет желаемого допуска.

tStart = 0          % Время началаКОНЕЦ = 5            % Время окончанияж = -у^2            % Производная y, поэтому y '= f (t, y (t)) = -y ^ 2                    % Решение этого ОДУ: y = 1 / (1 + t)y0 = 1              % Начальная позиция (т.е. y0 = y (tStart) = y (0) = 1)толерантность = 10^-11  % 10 желательна точностьmaxRows = 20                % Не позволяйте итерации продолжаться бесконечноinitialH = tStart - КОНЕЦ    % Выберите начальный размер шагаhaveWeFoundSolution = ложный % Удалось ли найти решение в пределах желаемого допуска? еще нет.час = initialH% Создайте 2D-матрицу размером maxRows на maxRows для хранения экстраполяций Ричардсона% Обратите внимание, что это будет нижняя треугольная матрица и что на самом деле не более двух строк% требуется в любой момент вычислений.А = zeroMatrix(maxRows, maxRows)% Вычислить верхний левый элемент матрицыА(1, 1) = Трапециевидный(ж, tStart, КОНЕЦ, час, y0)% Каждая строка матрицы требует одного вызова Trapezoidal% Этот цикл начинается с заполнения второй строки матрицы, поскольку первая строка была вычислена вышедля я = 1 : maxRows - 1 % Начиная с i = 1, повторять не более maxRows - 1 раз    час = час/2             % Половина предыдущего значения h, так как это начало новой строки        % Вызовите трапециевидную функцию с этим новым меньшим размером шага    А(я + 1, 1) = Трапециевидный(ж, tStart, КОНЕЦ, час, y0)    для j = 1: i% Пройдите по строке, пока не достигнете диагонали        % Используйте только что вычисленное значение (например, A (i + 1, j)) и элемент из        % строка над ним (т.е.A (i, j)) для вычисления следующей экстраполяции Ричардсона             А(я + 1, j + 1) = ((4^j).*А(я + 1, j) - А(я, j))/(4^j - 1);    конец% После выхода из указанного выше внутреннего цикла был вычислен диагональный элемент строки i + 1.    % Этот диагональный элемент представляет собой последнюю экстраполяцию Ричардсона, которая должна быть вычислена.    % Разница между этой экстраполяцией и последней экстраполяцией строки i является хорошей    % индикации ошибки.    если (абсолютная величина(А(я + 1, я + 1) - А(я, я)) < толерантность)  % Если результат в пределах допуска        Распечатать("у (5) =", А(я + 1, я + 1))                      % Показать результат экстраполяции Ричардсона        haveWeFoundSolution = правда        перерыв % Готово, выходите из цикла    конецконецесли (haveWeFoundSolution == ложный)   % Если бы мы не смогли найти решение в пределах желаемого допуска    Распечатать(«Предупреждение: невозможно найти решение в пределах желаемого допуска», толерантность);    Распечатать(«Последняя вычисленная экстраполяция была», А(maxRows, maxRows))конец

Смотрите также

использованная литература

^ Ричардсон, Л.Ф. (1911). «Приближенное арифметическое решение конечными разностями физических задач, включая дифференциальные уравнения, с приложением к напряжениям в каменной дамбе». Философские труды Королевского общества A. 210 (459–470): 307–357. Дои:10.1098 / рста.1911.0009.
^ Ричардсон, Л.Ф.; Гонт, Дж. А. (1927). «Отсроченный подход к пределу». Философские труды Королевского общества A. 226 (636–646): 299–349. Дои:10.1098 / рста.1927.0008.
^ Страница 126 из Биркофф, Гарретт; Джан-Карло Рота (1978). Обыкновенные дифференциальные уравнения (3-е изд.). Джон Уайли и сыновья. ISBN 0-471-07411-X. OCLC 4379402.

Методы экстраполяции. Теория и практика К. Брезински и М. Редиво Загля, Северная Голландия, 1991.
Иван Димов, Захари Златев, Иштван Фараго, Агнес Хаваси: Экстраполяция Ричардсона: практические аспекты и приложения, Walter de Gruyter GmbH & Co KG, ISBN 9783110533002 (2017).

внешние ссылки

Основные методы численной экстраполяции с приложениями, mit.edu
Ричардсон-Экстраполяция
Экстраполяция Ричардсона на сайте Роберта Исраэля (Университет Британской Колумбии)
Код Matlab для общей экстраполяции Ричардсона.
Код Юлии для общей экстраполяции Ричардсона.

[1] Ричардсон, Л.Ф. (1911). «Приближенное арифметическое решение конечными разностями физических задач, включая дифференциальные уравнения, с приложением к напряжениям в каменной дамбе». Философские труды Королевского общества A. 210 (459–470): 307–357. Дои:10.1098 / рста.1911.0009.

[2] Ричардсон, Л.Ф.; Гонт, Дж. А. (1927). «Отсроченный подход к пределу». Философские труды Королевского общества A. 226 (636–646): 299–349. Дои:10.1098 / рста.1927.0008.

[3] Страница 126 из Биркофф, Гарретт; Джан-Карло Рота (1978). Обыкновенные дифференциальные уравнения (3-е изд.). Джон Уайли и сыновья. ISBN 0-471-07411-X. OCLC 4379402.

[1]

[2]

[3]