метод Стеффенсена

В численном анализе метод Стеффенсена представляет собой итерационный метод поиска корней, названный в честь Йохана Фредерика Стеффенсена , который аналогичен методу Ньютона , но имеет определенные ситуационные преимущества. В частности, метод Стеффенсена достигает аналогичной квадратичной сходимости , но без использования производных , как это требуется для метода Ньютона .

Простое описание

Самая простая форма формулы метода Стеффенсена возникает, когда она используется для нахождения нуля действительной функции. $\ f\,;$ то есть найти реальную стоимость $x_{\star }$ это удовлетворяет $\ f(x_{\star })=0~.$ Рядом с решением $\ x_{\star }\,,$ функция $\ f\$ должно приблизительно удовлетворять $\ -1<f'(x_{\star })<0\,;$ это условие делает $\ f\$ адекватен в качестве корректирующей функции для $~x~$ для поиска собственного решения, хотя от него не требуется эффективной работы. Для некоторых функций метод Стеффенсена может работать, даже если это условие не выполняется, но в таком случае начальное значение $\ x_{0}\$ должно быть очень близко к реальному решению $\ x_{\star }\,,$ и сходимость к решению может быть медленной. Корректировка размера шага метода, упомянутая ниже, может улучшить сходимость в некоторых из этих случаев.

При адекватном стартовом значении $\ x_{0}\,,$ последовательность значений $\ x_{0},\ x_{1},\ x_{2},\dots ,\ x_{n},\ \dots \$ можно сгенерировать по приведенной ниже формуле. Когда это работает, каждое значение в последовательности намного ближе к решению. $\ x_{\star }\$ чем предыдущее значение. Значение $\ x_{n}\$ с текущего шага генерирует значение $\ x_{n+1}\$ для следующего шага по этой формуле:

$\ x_{n+1}=x_{n}-{\frac {\,f(x_{n})\,}{g(x_{n})}}\$

для $\ n=0,1,2,3,...\,;$ где функция наклона $\ g(x)\$ является композицией исходной функции $\ f\$ заданной по следующей формуле:

$\ g(x)={\frac {\,f{\bigl (}x+f(x){\bigr )}\,}{f(x)}}-1\$

или, возможно, более четко,

$\ g(x)~=~{\frac {\ f(x+h)-f(x)\ }{h}}\ ~~\approx ~~{\frac {\ \operatorname {d} f(x)\ }{\operatorname {d} x}}~\equiv ~f'(x)\ ,$

где $\ h=f(x)\$ — размер шага между последней точкой итерации, $\ x\ ,$ и вспомогательная точка, расположенная в $\ x+h~.$

Технически функция $\ g\$ порядка называется разделенной разностью первого $\ f\$ между этими двумя точками (это разделенная разность либо прямого , либо обратного типа, в зависимости от знака $\ h\$). Практически это усредненное значение наклона $\ f'\$ функции $\ f\$ между последней точкой последовательности $\ \left(x,y\right)={\bigl (}x_{n},\ f\left(x_{n}\right){\bigr )}\$ и вспомогательная точка в $\ {\bigl (}x,y{\bigr )}={\bigl (}\ x_{n}+h,\ f\left(x_{n}+h\right)\ {\bigr )}\ ,$ с размером шага (и направлением) $\ h=f(x_{n})~.$

Поскольку значение $\ g\$ является приближением для $\ f'\ ,$ его значение можно дополнительно проверить, чтобы увидеть, соответствует ли оно условию $\ -1<g<0\$ что требуется от $\ f'\ ,$ чтобы гарантировать сходимость алгоритма Стеффенсена. Хотя небольшое несоответствие не обязательно может быть ужасным, любое значительное отклонение от условия предупреждает, что метод Стеффенсена может потерпеть неудачу, и оправдано временное использование некоторого запасного алгоритма (например, более надежного алгоритма Иллинойса или простого regula falsi ).

Это только с целью найти $\ h\$ для этой вспомогательной точки значение функции $\ f\$ должна быть адекватной коррекцией, чтобы приблизиться к собственному решению и по этой причине выполнить требование, $\ -1<f'(x_{\star })<0~.$ Для всех остальных частей расчета метод Стеффенсена требует только функции $\ f\$ быть непрерывным и иметь близкое решение. Несколько скромных модификаций ступени $\ h\$ в формуле наклона $\ g\$ существуют, например, умножив его на ⁠ 1 / 2 ⁠ или ⁠ 3 / 4 ⁠ для размещения функций $\ f\$ что не совсем соответствует требованию.

Преимущества и недостатки

Основное преимущество метода Стеффенсена состоит в том, что он имеет квадратичную сходимость. как метод Ньютона , то есть оба метода находят корни уравнения $~f~$ так же «быстро». В данном случае быстро означает, что для обоих методов количество правильных цифр в ответе удваивается с каждым шагом. Но формула метода Ньютона требует вычисления производной функции. $~f'~$ а также функция $~f~,$ тогда как метод Стеффенсена требует только $~f~$ сам. Это важно, когда производный инструмент недоступен легко или эффективно.

Платой за быструю сходимость является двойная оценка функции: обе $~f(x_{n})~$ и $~f(x_{n}+h)~$ необходимо рассчитать, что может занять много времени, если $~f~$ это сложная функция. Для сравнения, метод секущего требует только одной оценки функции на шаг. Метод секущего увеличивает количество правильных цифр «всего» примерно в 1,6 раза за шаг, но за заданное время можно выполнить вдвое больше шагов метода секущего. Поскольку метод секущей может выполнять за то же время вдвое больше шагов, чем метод Стеффенсена, на практике метод секущего фактически сходится быстрее, чем метод Стеффенсена, когда оба алгоритма успешны: метод секущего достигает коэффициента примерно (1,6). ² ≈ в 2,6 раза больше цифр на каждые два шага (две оценки функции) по сравнению с коэффициентом Стеффенсена, равным 2 на каждый шаг (две оценки функции).

Как и в большинстве других итеративных алгоритмов поиска корня , решающей слабостью метода Стеффенсена является выбор «адекватного» начального значения. $~x_{0}~.$ Если значение $~x_{0}~$ не «достаточно близко» к фактическому решению $~x_{\star }~,$ метод может дать сбой и последовательность значений $~x_{0},\,x_{1},\,x_{2},\,x_{3},\,\dots ~$ может либо хаотично переключаться между двумя крайностями, либо расходиться до бесконечности, либо и то, и другое.

Вывод с использованием процесса Эйткена дельта-квадрат.

Версию метода Стеффенсена, реализованную в коде MATLAB , показанном ниже, можно найти, используя процесс дельта-квадрата Эйткена для ускорения сходимости последовательности. Чтобы сравнить следующие формулы с формулами из раздела выше, обратите внимание, что $x_{n}=p\,-\,p_{n}~.$ Этот метод предполагает начинать с линейно сходящейся последовательности и увеличивает скорость сходимости этой последовательности. Если признаки $~p_{n},\,p_{n+1},\,p_{n+2}~$ согласен и $~p_{n}~$ «достаточно близко» к желаемому пределу последовательности $~p~,$ мы можем предположить следующее:

${\frac {\,p_{n+1}-p\,}{\,p_{n}-p\,}}~\approx ~{\frac {\,p_{n+2}-p\,}{\,p_{n+1}-p\,}}$

затем

$(p_{n+1}-p)^{2}~\approx ~(\,p_{n+2}-p\,)\,(\,p_{n}-p\,)~$

так

$p_{n+1}^{2}-2\,p_{n+1}\,p+p^{2}~\approx ~p_{n+2}\;p_{n}-(\,p_{n}+p_{n+2}\,)\,p+p^{2}$

и, следовательно,

$(\,p_{n+2}-2\,p_{n+1}+p_{n}\,)\,p~\approx ~p_{n+2}\,p_{n}-p_{n+1}^{2}~.$

Решение для желаемого предела последовательности $~p~$ дает:

$p~\approx ~{\frac {\,p_{n+2}\,p_{n}-p_{n+1}^{2}\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}~=~{\frac {\,p_{n}^{2}+p_{n}\,p_{n+2}+2\,p_{n}\,p_{n+1}-2\,p_{n}\,p_{n+1}-p_{n}^{2}-p_{n+1}^{2}\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}$

$=~{\frac {\,(\,p_{n}^{2}+p_{n}\,p_{n+2}-2\,p_{n}\,p_{n+1}\,)-(\,p_{n}^{2}-2\,p_{n}\,p_{n+1}+p_{n+1}^{2}\,)\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}$

$=~p_{n}-{\frac {\,(\,p_{n+1}-p_{n})^{2}\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}~,$

что приводит к более быстро сходящейся последовательности:

$p~\approx ~p_{n+3}~=~p_{n}-{\frac {\,(\,p_{n+1}-p_{n}\,)^{2}\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}~.$

Пример кода

В Матлабе

Вот источник реализации метода Стеффенсена в MATLAB .

function Steffensen(f, p0, tol)
% This function takes as inputs: a fixed point iteration function, f, 
% and initial guess to the fixed point, p0, and a tolerance, tol.
% The fixed point iteration function is assumed to be input as an
% inline function. 
% This function will calculate and return the fixed point, p, 
% that makes the expression f(x) = p true to within the desired 
% tolerance, tol.

format compact   % This shortens the output.
format long      % This prints more decimal places.

for i = 1:1000   % get ready to do a large, but finite, number of iterations.
                 % This is so that if the method fails to converge, we won't
                 % be stuck in an infinite loop.
    p1 = f(p0);  % calculate the next two guesses for the fixed point.
    p2 = f(p1);
    p = p0-(p1-p0)^2/(p2-2*p1+p0) % use Aitken's delta squared method to
                                  % find a better approximation to p0.
    if abs(p - p0) < tol  % test to see if we are within tolerance.
        break             % if we are, stop the iterations, we have our answer.
    end
    p0 = p;               % update p0 for the next iteration.
end
if abs(p - p0) > tol      % If we fail to meet the tolerance, we output a
                          % message of failure.
    'failed to converge in 1000 iterations.'
end

На Python

Вот исходный код реализации метода Стеффенсена на Python .

from typing import Callable, Iterator
Func = Callable[[float], float]

def g(f: Func, x: float, fx: float) -> Func:
    """First-order divided difference function.

    Arguments:
        f: Function input to g
        x: Point at which to evaluate g
        fx: Function f evaluated at x 
    """
    return lambda x: f(x + fx) / fx - 1

def steff(f: Func, x: float) -> Iterator[float]:
    """Steffenson algorithm for finding roots.

    This recursive generator yields the x_{n+1} value first then, when the generator iterates,
    it yields x_{n+2} from the next level of recursion.

    Arguments:
        f: Function whose root we are searching for
        x: Starting value upon first call, each level n that the function recurses x is x_n
    """
    while True:    
        fx = f(x)
        gx = g(f, x, fx)(x)
        if gx == 0:
            break
        else:
            x = x - fx / gx    # Update to x_{n+1}
            yield x            # Yield value

Обобщение на банахово пространство

Метод Стеффенсена также можно использовать для поиска входных данных. $\ x=x_{\star }\$ для другого типа функции $\ F\$ который производит вывод такой же, как и ввод: $\ x_{\star }=F(x_{\star })\$ за особую ценность $\ x_{\star }~.$ Такие решения, как $\ x_{\star }\$ называются фиксированными точками . Многие из этих функций можно использовать для поиска собственных решений путем многократного повторного использования результата обратно в качестве входных данных, но скорость сходимости может быть медленной или функция может вообще не сходиться, в зависимости от отдельной функции. Метод Стеффенсена ускоряет эту сходимость, делая ее квадратичной .

Для ориентации корневая функция $\ f\$ а функции фиксированной точки просто связаны соотношением $\ F(x)=x+\varepsilon \ f(x)\ ,$ где $\ \varepsilon \$ является некоторой скалярной константой, достаточно малой по величине, чтобы сделать $\ F\$ стабилен при итерации, но достаточно велик для нелинейности функции $\ f\$ быть заметным. Все вопросы, связанные с более общим банаховым пространством и базовыми действительными числами, на мгновение игнорируются ради сравнения.

Этот метод нахождения неподвижных точек вещественной функции обобщен на функции $\ F:X\to X\$ которые отображают банахово пространство $\ X\$ на себя или даже в более общем плане $\ F:X\to Y\$ это отображение из одного банахова пространства $\ X\$ в другое банахово пространство $\ Y~.$ Обобщенный метод предполагает, что ограниченных линейных семейство операторов $\ \{\;G(u,v):u,v\in X\;\}\$ связанный с $\ u\$ и $\ v\$ можно придумать, что (локально) удовлетворяет этому условию:

$\ F\left(u\right)-F\left(v\right)=G\left(u,v\right)\,{\bigl [}\,u-v\,{\bigr ]}\ \qquad \qquad \qquad \qquad$

уравнение ( ⸎ )

Если деление возможно в банаховом пространстве , то линейный оператор $\ G\$ можно получить из

$\ G\left(u,v\right)={\bigl [}\ F\left(u\right)-F\left(v\right)\ {\bigr ]}{\bigl [}\ u-v\ {\bigr ]}^{-1}\ ,$

что может дать некоторое представление: выраженный таким образом линейный оператор $\ G\$ легче рассматривать как сложную версию разделенной разности. $\ g\$ обсуждалось в первом разделе выше. Форма частного показана здесь только для ориентации; это не требуется само по себе . Обратите также внимание, что разделение внутри банахового пространства не является необходимым для жизнеспособности разработанного метода Стеффенсена; единственное требование состоит в том, чтобы оператор $\ G\$ удовлетворяют отмеченному короной уравнению , , ( ⸎ ) .

Для основной функции вещественного числа $\ f\ ,$ Как указано в первом разделе, функция просто принимает и выдает действительные числа. Там функция $\ g\$ это разделенная разница . В обобщенном виде здесь оператор $\ G\$ является аналогом разделенной разности для использования в банаховом пространстве . Оператор $\ G\$ примерно эквивалентно матрице , все элементы которой являются функциями векторных аргументов $\ u\$ и $\ v~.$

В этом случае метод Стеффенсена очень похож на метод Ньютона, за исключением того, что он использует разделенную разность. $\ G{\bigl (}\,F\left(x\right),\,x\,{\bigr )}\$ вместо производной $\ F'(x)~.$ Обратите внимание, что для аргументов $\ x\$ близко к некоторой фиксированной точке $\ x_{\star }\ ,$ функции с фиксированной точкой $\ F\$ и их линейные операторы $\ G\$ соответствие условию, отмеченному ( ⸎ ) , $\ F'(x)\approx G{\bigl (}\,F\left(x\right),\,x\,{\bigr )}\approx I\$ где $\ I\$ является идентификационным оператором.

В случае, когда деление возможно в банаховом пространстве, обобщенная итерационная формула имеет вид

$x_{n+1}=x_{n}+{\Bigl [}I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}{\Bigr ]}^{-1}{\Bigl [}F\left(x_{n}\right)-x_{n}{\Bigr ]}\ ,$

для $\ n=1,\,2,\,3,\,...~.$ В более общем случае, когда деление может оказаться невозможным, итерационная формула требует поиска решения $\ x_{n+1}\$ близко к $\ x_{n}\$ для чего

${\Bigl [}I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}{\Bigr ]}{\Bigl [}x_{n+1}-x_{n}{\Bigr ]}=F\left(x_{n}\right)-x_{n}~.$

Аналогично можно искать решение $\ x_{n+1}\$ в несколько уменьшенную форму

${\Bigl [}I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}{\Bigr ]}x_{n+1}={\Bigl [}F\left(x_{n}\right)-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ x_{n}{\Bigr ]}\ ,$

при этом все значения внутри квадратных скобок не зависят от $\ x_{n+1}\ :$ Члены в скобках зависят только от $\ x_{n}~.$ Однако имейте в виду, что вторая форма может быть не такой численно стабильной, как первая: поскольку первая форма включает в себя поиск значения для (надеюсь) небольшой разницы, в числовом отношении более вероятно избежать чрезмерно больших или беспорядочных изменений в итерируемом методе. ценить $\ x_{n}~.$

Если линейный оператор $~G~$ удовлетворяет

${\Bigl \|}G\left(u,v\right)-G\left(x,y\right){\Bigr \|}\leq k{\biggl (}{\Bigl \|}u-x{\Bigr \|}+{\Bigr \|}v-y{\Bigr \|}{\biggr )}~$

для некоторой положительной действительной константы $\ k~,$ то метод сходится квадратично к фиксированной точке $\ F\$ если начальное приближение $\ x_{0}\$ «достаточно близко» к желаемому решению $\ x_{\star }\ ,$ это удовлетворяет $\ x_{\star }=F(x_{\star })~.$