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

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

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

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

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

Обозначения

Позволять $A_{0}(h)$ быть приближением $A^{*}$(точное значение), которое зависит от размера шага h (где ${\textstyle 0<h<1}$) с формулой ошибки вида

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

где $a_{i}$ неизвестные константы, а $k_{i}$ являются известными константами такими, что $h^{k_{i}}>h^{k_{i+1}}$. Более того, $O(h^{k_{i}})$ представляет собой усечения ошибку $A_{i}(h)$ приближение такое, что $A^{*}=A_{i}(h)+O(h^{k_{i}}).$ Аналогично, в $A^{*}=A_{i}(h)+O(h^{k_{i}}),$ приближение $A_{i}(h)$ Говорят, что это $O(h^{k_{i}})$ приближение.

Обратите внимание, что при упрощении с помощью обозначения Big O следующие формулы эквивалентны:

$${\begin{aligned}A^{*}&=A_{0}(h)+a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+\cdots \\A^{*}&=A_{0}(h)+a_{0}h^{k_{0}}+O(h^{k_{1}})\\A^{*}&=A_{0}(h)+O(h^{k_{0}})\end{aligned}}$$

Цель

Экстраполяция Ричардсона – это процесс, позволяющий найти лучшее приближение $A^{*}$ изменив формулу ошибки с $A^{*}=A_{0}(h)+O(h^{k_{0}})$ к $A^{*}=A_{1}(h)+O(h^{k_{1}}).$ Поэтому, заменив $A_{0}(h)$ с $A_{1}(h)$ ошибка усечения уменьшилась с $O(h^{k_{0}})$ к $O(h^{k_{1}})$ для того же размера шага $h$. Происходит общая закономерность, при которой $A_{i}(h)$ является более точной оценкой, чем $A_{j}(h)$ когда $i>j$. Благодаря этому процессу мы достигли лучшего приближения $A^{*}$ вычитая самый большой член ошибки, которая была $O(h^{k_{0}})$. Этот процесс можно повторить, чтобы удалить больше ошибок и получить еще лучшие приближения.

Процесс

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

$$A^{*}=A_{0}(h)+a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+O(h^{k_{3}})$$

( 1 )

$$A^{*}=A_{0}\!\left({\frac {h}{t}}\right)+a_{0}\left({\frac {h}{t}}\right)^{k_{0}}+a_{1}\left({\frac {h}{t}}\right)^{k_{1}}+a_{2}\left({\frac {h}{t}}\right)^{k_{2}}+O(h^{k_{3}})$$

( 2 )

Чтобы улучшить наше приближение из $O(h^{k_{0}})$ к $O(h^{k_{1}})$ удалив первый член ошибки, мы умножаем уравнение 2 на $t^{k_{0}}$ и вычтем уравнение 1, чтобы получить

$$(t^{k_{0}}-1)A^{*}={\bigg [}t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h){\bigg ]}+{\bigg (}t^{k_{0}}a_{1}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{1}}-a_{1}h^{k_{1}}{\bigg )}+{\bigg (}t^{k_{0}}a_{2}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{2}}-a_{2}h^{k_{2}}{\bigg )}+O(h^{k_{3}}).$$

Это умножение и вычитание выполнялось потому, что ${\textstyle {\big [}t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h){\big ]}}$ это $O(h^{k_{1}})$ приближение $(t^{k_{0}}-1)A^{*}$. Мы можем решить нашу текущую формулу для $A^{*}$ дать

$$A^{*}={\frac {{\bigg [}t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h){\bigg ]}}{t^{k_{0}}-1}}+{\frac {{\bigg (}t^{k_{0}}a_{1}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{1}}-a_{1}h^{k_{1}}{\bigg )}}{t^{k_{0}}-1}}+{\frac {{\bigg (}t^{k_{0}}a_{2}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{2}}-a_{2}h^{k_{2}}{\bigg )}}{t^{k_{0}}-1}}+O(h^{k_{3}})$$

который можно записать как $A^{*}=A_{1}(h)+O(h^{k_{1}})$ установив

$$A_{1}(h)={\frac {t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h)}{t^{k_{0}}-1}}.$$

Рекуррентное отношение

Общее рекуррентное соотношение можно определить для аппроксимаций следующим образом:

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

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

$$A^{*}=A_{i+1}(h)+O(h^{k_{i+1}}).$$

Свойства

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

Кроме того, общую формулу можно использовать для оценки $k_{0}$ (поведение размера шага ведущего порядка при ошибке усечения ), когда ни его значение, ни $A^{*}$ известно априори . Такой метод может быть полезен для количественного определения неизвестной скорости сходимости . Учитывая приближения $A^{*}$ из трех различных размеров шага $h$, $h/t$, и $h/s$, точное соотношение

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

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

$$A_{i}\left({\frac {h}{t}}\right)+{\frac {A_{i}\left({\frac {h}{t}}\right)-A_{i}(h)}{t^{k_{0}}-1}}\approx A_{i}\left({\frac {h}{s}}\right)+{\frac {A_{i}\left({\frac {h}{s}}\right)-A_{i}(h)}{s^{k_{0}}-1}}$$

которое можно решить численно, чтобы оценить $k_{0}$ для некоторого произвольного допустимого выбора $h$, $s$, и $t$.

Как $t\neq 1$, если $t>0$ и $s$ выбирается так, что $s=t^{2}$, это приближенное соотношение сводится к квадратному уравнению относительно $t^{k_{0}}$, что легко решается для $k_{0}$ с точки зрения $h$ и $t$.

Пример

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

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

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

$$R(h,t):={\frac {t^{n}A(h/t)-A(h)}{t^{n}-1}}$$

где $h$ и ${\frac {h}{t}}$ это два разных размера шага.

Затем

$$R(h,t)={\frac {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}).$$

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

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

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

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

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

tStart = 0          % Starting time
tEnd = 5            % Ending time
f = -y^2            % The derivative of y, so y' = f(t, y(t)) = -y^2
                    % The solution to this ODE is y = 1/(1 + t)
y0 = 1              % The initial position (i.e. y0 = y(tStart) = y(0) = 1)
tolerance = 10^-11  % 10 digit accuracy is desired

% Don't allow the iteration to continue indefinitely
maxRows = 20
% Pick an initial step size
initialH = tStart - tEnd
% Were we able to find the solution to within the desired tolerance? not yet.
haveWeFoundSolution = false

h = initialH

% Create a 2D matrix of size maxRows by maxRows to hold the Richardson extrapolates
% Note that this will be a lower triangular matrix and that at most two rows are actually
% needed at any time in the computation.
A = zeroMatrix(maxRows, maxRows)

% Compute the top left element of the matrix.
% The first row of this (lower triangular) matrix has now been filled.
A(1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

% Each row of the matrix requires one call to Trapezoidal
% This loops starts by filling the second row of the matrix,
% since the first row was computed above
for i = 1 : maxRows - 1 % Starting at i = 1, iterate at most maxRows - 1 times
    % Halve the previous value of h since this is the start of a new row.
    h = h/2

    % Starting filling row i+1 from the left by calling
    % the Trapezoidal function with this new smaller step size
    A(i + 1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

    % Go across this current (i+1)-th row until the diagonal is reached
    for j = 1 : i
        % To compute A(i + 1, j + 1), which is the next Richardson extrapolate, 
        % use the most recently computed value (i.e. A(i + 1, j))
        % and the value from the row above it (i.e. A(i, j)).

        A(i + 1, j + 1) = ((4^j).*A(i + 1, j) - A(i, j))/(4^j - 1);
    end
    
    % After leaving the above inner loop, the diagonal element of row i + 1 has been computed
    % This diagonal element is the latest Richardson extrapolate to be computed.
    % The difference between this extrapolate and the last extrapolate of row i is a good
    % indication of the error.
    if (absoluteValue(A(i + 1, i + 1) - A(i, i)) < tolerance)  % If the result is within tolerance
        % Display the result of the Richardson extrapolation
        print("y = ", A(i + 1, i + 1))
        haveWeFoundSolution = true
        % Done, so leave the loop
        break
    end
end

% If we were not able to find a solution to within the desired tolerance
if (not haveWeFoundSolution)
    print("Warning: Not able to find solution to within the desired tolerance of ", tolerance);
    print("The last computed extrapolate was ", A(maxRows, maxRows))
end

См. также

• Дельта-квадратичный процесс Эйткена
• Такэбе Кенко
• Итерация Ричардсона

Ссылки

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

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

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