Экстраполяция Ричардсона - Richardson extrapolation
В численном анализе , Экстраполяция Ричардсона является ускорение последовательности метод , используемая для улучшения скорости сходимости в виде последовательности оценок некоторого значения . По сути, учитывая значение для нескольких значений , мы можем оценить , экстраполируя оценки на . Он назван в честь Льюиса Фрая Ричардсона , который представил эту технику в начале 20-го века, хотя идея была уже известна Христиану Гюйгенсу в его вычислении π . По словам Биркгофа и Роты , «его полезность для практических вычислений трудно переоценить».
Практические применения экстраполяции Ричардсона включают интегрирование Ромберга , которое применяет экстраполяцию Ричардсона к правилу трапеций , и алгоритм Булирша – Стоера для решения обыкновенных дифференциальных уравнений.
Пример экстраполяции Ричардсона
Предположим, что мы хотим приблизиться , и у нас есть метод, который зависит от небольшого параметра таким образом, что
Определим новую функцию
где и - два различных размера шага.
потом
называется Ричардсон экстраполяцией из А ( ч ), и имеет оценку погрешности более высокого порядка по сравнению с .
Очень часто гораздо проще получить заданную точность, используя R (h), а не A (h ') с гораздо меньшим h' . Где A (h ') может вызвать проблемы из-за ограниченной точности ( ошибки округления ) и / или из-за увеличения количества необходимых вычислений (см. Примеры ниже).
Общая формула
Позвольте быть приближением (точное значение), которое зависит от положительного размера шага h с формулой ошибки вида
где a i - неизвестные константы, а k i - известные константы, такие что h k i > h k i + 1 .
Кроме того, k 0 - это поведение размера шага ведущего порядка ошибки усечения, как
Искомое точное значение может быть выражено как
который можно упростить с помощью нотации Big O до
Используя размеры шага h и некоторую константу t , две формулы для A :
Умножение второго уравнения на t k 0 и вычитание первого уравнения дает
который можно решить, чтобы дать
Таким образом, ошибка усечения была уменьшена до . Это в отличие от , где ошибка усечения находится за тот же размер шага
Благодаря этому процессу мы достигли лучшего приближения к A , вычтя наибольший член ошибки, который был O ( h k 0 ). Этот процесс можно повторить, чтобы удалить больше ошибок, чтобы получить еще лучшие приближения.
Общее рекуррентное соотношение, начинающееся с, может быть определено для приближений формулой
где удовлетворяет
- .
Экстраполяцию Ричардсона можно рассматривать как преобразование линейной последовательности .
Кроме того, общая формула может использоваться для оценки k 0 (поведение размера шага ведущего порядка ошибки усечения ), когда ни его значение, ни A * (точное значение) не известны априори . Такой метод может быть полезен для количественной оценки неизвестной скорости сходимости . Учитывая приближения A из трех различных размеров шага h , h / t и h / s , точное соотношение
дает приблизительное соотношение (обратите внимание, что обозначения здесь могут вызвать небольшую путаницу, два O, появляющиеся в приведенном выше уравнении, указывают только на поведение размера шага ведущего порядка, но их явные формы различны, и, следовательно, сокращение двух членов O является приблизительно действительный)
которое можно решить численно, чтобы оценить k 0 для некоторых произвольных выборов h , s и t .
Пример кода псевдокода для экстраполяции Ричардсона
Следующий псевдокод в стиле MATLAB демонстрирует Экстраполяция Ричардсона , чтобы помочь решить ОДУ , с методом трапеций . В этом примере мы уменьшаем вдвое размер шага на каждой итерации, и поэтому в приведенном выше обсуждении это будет . Ошибка трапециевидного метода может быть выражена в терминах нечетных степеней, так что ошибку на нескольких шагах можно выразить в четных степенях; это приводит нас к возведению во вторую степень и к степеням псевдокода. Мы хотим найти значение , которое имеет точное решение, поскольку точное решение ОДУ равно . Этот псевдокод предполагает, что существует вызываемая функция, которая пытается вычислить , выполняя трапециевидный метод для функции с начальной точкой и размером шага .
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
maxRows = 20 % Don't allow the iteration to continue indefinitely
initialH = tStart - tEnd % Pick an initial step size
haveWeFoundSolution = false % Were we able to find the solution to within the desired tolerance? not yet.
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
h = h/2 % Halve the previous value of h since this is the start of a new row.
% 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)
for j = 1 : i % Go across this current (i+1)-th row until the diagonal is reached
% 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
print("y = ", A(i + 1, i + 1)) % Display the result of the Richardson extrapolation
haveWeFoundSolution = true
break % Done, so leave the loop
end
end
if (haveWeFoundSolution == false) % If we were not able to find a solution to within the desired tolerance
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 для обобщенной экстраполяции Ричардсона.
- Код Джулии для общей экстраполяции Ричардсона.