Метод Ньютона

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

Описание методаПравить

ОбоснованиеПравить

Чтобы численно решить уравнение $f(x)=0$$f(x)=0$  методом простой итерации, его необходимо привести к эквивалентному уравнению: $x=\mathrm{\phi <!--\; \phi \; -->}(x)$$x=\varphi(x)$ , где $\mathrm{\phi <!--\; \phi \; -->}$$\varphi$  — сжимающее отображение.

Для наилучшей сходимости метода в точке очередного приближения $x^{\ast <!--\; \ast \; -->}$$x^{*}$  должно выполняться условие ${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x^{\ast <!--\; \ast \; -->})=0$$\varphi '(x^{*})=0$ . Решение данного уравнения ищут в виде $\mathrm{\phi <!--\; \phi \; -->}(x)=x+\mathrm{\alpha <!--\; \alpha \; -->}(x)f(x)$$\varphi (x)=x+\alpha (x)f(x)$ , тогда:

${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x^{\ast <!--\; \ast \; -->})=1+{\mathrm{\alpha <!--\; \alpha \; -->}}^{\prime }(x^{\ast <!--\; \ast \; -->})f(x^{\ast <!--\; \ast \; -->})+\mathrm{\alpha <!--\; \alpha \; -->}(x^{\ast <!--\; \ast \; -->})f^{\prime }(x^{\ast <!--\; \ast \; -->})=0.$$\varphi '(x^{*})=1+\alpha '(x^{*})f(x^{*})+\alpha (x^{*})f'(x^{*})=0.$ 

В предположении, что точка приближения «достаточно близка» к корню $\stackrel{~<!--\; ~\; -->}{x}$${\tilde {x}}$  и что заданная функция непрерывна $(f(x^{\ast <!--\; \ast \; -->})\approx <!--\; \approx \; -->f(\stackrel{~<!--\; ~\; -->}{x})=0)$$(f(x^{*})\approx f({\tilde {x}})=0)$ , окончательная формула для $\mathrm{\alpha <!--\; \alpha \; -->}(x)$$\alpha(x)$  такова:

$\mathrm{\alpha <!--\; \alpha \; -->}(x)=-<!--\; -\; -->\frac{1}{f^{\prime }(x)}.$$\alpha (x)=-{\frac {1}{f'(x)}}.$ 

С учётом этого функция $\mathrm{\phi <!--\; \phi \; -->}(x)$$\varphi (x)$  определяется:

$\mathrm{\phi <!--\; \phi \; -->}(x)=x-<!--\; -\; -->\frac{f(x)}{f^{\prime }(x)}.$$\varphi (x)=x-{\frac {f(x)}{f'(x)}}.$ 

При некоторых условиях эта функция в окрестности корня осуществляет сжимающее отображение.

Доказательство

Пусть дана функция вещественного переменного дважды непрерывно дифференцируемая в своей области определения, производная которой нигде не обращается в нуль:

$f(x):<!--\; :\; -->\mathbb{X}\to <!--\; \to \; -->\mathbb{R},f(x)\in <!--\; \in \; -->{\mathrm{C}}^2(\mathbb{X});\mathrm{\forall <!--\; \forall \; -->}x\in <!--\; \in \; -->\mathbb{X}{f}^{\prime }(x)\ne <!--\; \ne \; -->0.$$\scriptstyle {f(x)\colon {\mathbb {X}}\to \mathbb{R} ,\;f(x)\in {\mathrm {C}}^{2}({\mathbb {X}});\quad \forall x\in {\mathbb {X}}\;f'(x)\neq 0.}$ 

И необходимо доказать, что функция $\mathrm{\phi <!--\; \phi \; -->}(x)=x-<!--\; -\; -->\frac{f(x)}{f^{\prime }(x)}$$\scriptstyle {\varphi (x)=x-{\frac {f(x)}{f'(x)}}}$  осуществляет сжимающее отображение вблизи корня уравнения $f(x)=0$$\scriptstyle {f(x)=0}$ .

В силу непрерывной дифференцируемости функции $f(x)$$\scriptstyle {f(x)}$  и неравенства нулю её первой производной $\mathrm{\phi <!--\; \phi \; -->}(x)$$\scriptstyle {\varphi (x)}$  непрерывна.

Производная ${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x)$$\scriptstyle {\varphi '(x)}$  равна:

${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x)=\frac{f(x)f^{″}(x)}{{\left(f^{\prime }(x)\right)}^2}.$$\scriptstyle {\varphi '(x)={\frac {f(x)f''(x)}{\left(f'(x)\right)^{2}}}.}$ 

В условиях, наложенных на $f(x)$$\scriptstyle {f(x)}$ , она также непрерывна. Пусть $\stackrel{~<!--\; ~\; -->}{x}$$\scriptstyle {{\tilde {x}}}$  — искомый корень уравнения: $f(\stackrel{~<!--\; ~\; -->}{x})=0$$\scriptstyle {f({\tilde {x}})=0}$ , следовательно в его окрестности ${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x)\approx <!--\; \approx \; -->0$$\scriptstyle {\varphi '(x)\approx 0}$ :

$\scriptstyle {\forall \varepsilon \colon 0<\varepsilon <1,\;\exists \delta >0\;\forall x\in {\mathbb {X}}\;|x-{\tilde {x}}|<\delta \colon |\varphi '(x)-0|<\varepsilon .}$ 

Тогда согласно теореме Лагранжа:

$\scriptstyle {\forall x_{1},\;x_{2}\in {\mathrm {U}}_{\delta }({\tilde {x}})\;\exists \xi \in {\mathrm {U}}_{\delta }({\tilde {x}})\colon |\varphi (x_{1})-\varphi (x_{2})|=|\varphi '(\xi )||x_{1}-x_{2}|<\varepsilon |x_{1}-x_{2}|.}$ 

В силу того, что $\mathrm{\phi <!--\; \phi \; -->}(\stackrel{~<!--\; ~\; -->}{x})=\stackrel{~<!--\; ~\; -->}{x}$$\scriptstyle {\varphi ({\tilde {x}})={\tilde {x}}}$  в этой же дельта окрестности выполняется:

$\scriptstyle {\forall x\in U_{{\delta }}({\tilde {x}})\colon \;|\varphi (x)-{\tilde {x}}|<\varepsilon |x-{\tilde {x}}|.}$ 

Таким образом полученная функция $\mathrm{\phi <!--\; \phi \; -->}(x)$$\scriptstyle {\varphi (x)}$  в окрестности корня $U_{\mathrm{\delta <!--\; \delta \; -->}}(\stackrel{~<!--\; ~\; -->}{x})$$\scriptstyle {U_{\delta }({\tilde {x}})}$  осуществляет сжимающее отображение.■

В этом случае алгоритм нахождения численного решения уравнения $f(x)=0$$f(x)=0$  сводится к итерационной процедуре вычисления:

$x_{n+1}=x_n-<!--\; -\; -->\frac{f(x_n)}{f^{\prime }(x_n)}.$$x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}.$ 

По теореме Банаха последовательность приближений стремится к корню уравнения $f(x)=0$$f(x)=0$ .

 

Иллюстрация метода Ньютона (синим изображена функция $f(x)$$\scriptstyle {f(x)}$ , ноль которой необходимо найти, красным — касательная в точке очередного приближения $x_n$$\scriptstyle {x_{n}}$ ). Здесь мы можем увидеть, что последующее приближение $x_{n+1}$$\scriptstyle {x_{{n+1}}}$  лучше предыдущего $x_n$$\scriptstyle {x_{n}}$ .

Геометрическая интерпретацияПравить

Основная идея метода заключается в следующем: задаётся начальное приближение вблизи предположительного корня, после чего строится касательная к графику исследуемой функции в точке приближения, для которой находится пересечение с осью абсцисс. Эта точка берётся в качестве следующего приближения. И так далее, пока не будет достигнута необходимая точность.

Пусть 1) вещественнозначная функция $f(x):<!--\; :\; -->(a,b)\to <!--\; \to \; -->\mathbb{R}$$f(x)\colon (a,\,b)\to \mathbb {R}$  непрерывно дифференцируема на интервале $(a,b)$$(a,\,b)$  ;
2) существует искомая точка $x^{\ast <!--\; \ast \; -->}\in <!--\; \in \; -->(a,b)$$x^{*}\in (a,\,b)$  : $f(x^{\ast <!--\; \ast \; -->})=0$$f(x^{*})=0$  ;
3) существуют $C>0$$C > 0$  и $\mathrm{\delta <!--\; \delta \; -->}>0$$\delta >0$  такие, что
$|f^{\prime }(x)|⩾<!--\; ⩾\; -->C$$\vert f'(x)\vert \geqslant C$  для $x\in <!--\; \in \; -->(a,x^{\ast <!--\; \ast \; -->}-<!--\; -\; -->\mathrm{\delta <!--\; \delta \; -->}]\cup <!--\; \cup \; -->[x^{\ast <!--\; \ast \; -->}+\mathrm{\delta <!--\; \delta \; -->},b)$$x\in (a,\,x^{*}-\delta ]\cup [x^{*}+\delta ,\,b)$ 
и $f^{\prime }(x)\ne <!--\; \ne \; -->0$$f'(x)\neq 0$  для $x\in <!--\; \in \; -->(x^{\ast <!--\; \ast \; -->}-<!--\; -\; -->\mathrm{\delta <!--\; \delta \; -->},x^{\ast <!--\; \ast \; -->})\cup <!--\; \cup \; -->(x^{\ast <!--\; \ast \; -->},x^{\ast <!--\; \ast \; -->}+\mathrm{\delta <!--\; \delta \; -->})$$x\in (x^{*}-\delta ,\,x^{*})\cup (x^{*},\,x^{*}+\delta )$  ;
4) точка $x_n\in <!--\; \in \; -->(a,b)$$x_{n}\in (a,\,b)$  такова, что $f(x_n)\ne <!--\; \ne \; -->0$$f(x_{n})\neq 0$  .
Тогда формула итеративного приближения $x_n$$x_{n}$  к $x^{\ast <!--\; \ast \; -->}$$x^{{*}}$  может быть выведена из геометрического смысла касательной следующим образом:

$f^{\prime }(x_n)=\mathrm{t}\mathrm{g}{\mathrm{\alpha <!--\; \alpha \; -->}}_n=\frac{\mathrm{\Delta <!--\; \Delta \; -->}y}{\mathrm{\Delta <!--\; \Delta \; -->}x}=\frac{f(x_n)-<!--\; -\; -->0}{x_n-<!--\; -\; -->x_{n+1}}=\frac{0-<!--\; -\; -->f(x_n)}{x_{n+1}-<!--\; -\; -->x_n},$$f'(x_{n})=\mathrm {tg} \,\alpha _{n}={\frac {\Delta y}{\Delta x}}={\frac {f(x_{n})-0}{x_{n}-x_{n+1}}}={\frac {0-f(x_{n})}{x_{n+1}-x_{n}}},$ 

где ${\mathrm{\alpha <!--\; \alpha \; -->}}_n$$\alpha _{n}$  — угол наклона касательной прямой $y(x)=f(x_n)+(x-<!--\; -\; -->x_n)\cdot <!--\; \cdot \; -->\mathrm{t}\mathrm{g}{\mathrm{\alpha <!--\; \alpha \; -->}}_n$$y(x)=f(x_{n})+(x-x_{n})\cdot \mathrm {tg} \,\alpha _{n}$  к графику $f$$f$  в точке $(x_n;f(x_n))$$(x_{n};f(x_{n}))$  .

Следовательно (в уравнении касательной прямой полагаем $y(x_{n+1})=0$$y(x_{n+1})=0$ ) искомое выражение для $x_{n+1}$$x_{{n+1}}$  имеет вид :

$x_{n+1}=x_n-<!--\; -\; -->\frac{f(x_n)}{f^{\prime }(x_n)}.$$x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}.$ 

Если $x_{n+1}\in <!--\; \in \; -->(a,b)$$x_{n+1}\in (a,\,b)$ , то это значение можно использовать в качестве следующего приближения к $x^{\ast <!--\; \ast \; -->}$$x^{{*}}$  .

Если $x_{n+1}\notin <!--\; \notin \; -->(a,b)$$x_{n+1}\notin (a,\,b)$ , то имеет место «перелёт» (корень $x^{\ast <!--\; \ast \; -->}$$x^{{*}}$  лежит рядом с границей $(a,b)$$(a,\,b)$ ). В этом случае надо (воспользовавшись идеей метода половинного деления) заменять $x_{n+1}$$x_{{n+1}}$  на $\frac{x_n+x_{n+1}}{2}$${\frac {x_{n}+x_{n+1}}{2}}$  до тех пор, пока точка «не вернётся» в область поиска $(a,b)$$(a,\,b)$  .

Замечания. 1) Наличие непрерывной производной даёт возможность строить непрерывно меняющуюся касательную на всей области поиска решения $(a,b)$$(a,\,b)\;$  .
2) Случаи граничного (в точке $a$$a$  или в точке $b$$b$ ) расположения искомого решения $x^{\ast <!--\; \ast \; -->}$$x^{{*}}$  рассматриваются аналогичным образом.
3) С геометрической точки зрения равенство $f^{\prime }(x_n)=0$$f'(x_{n})=0$  означает, что касательная прямая к графику $f$$f$  в точке $(x_n;f(x_n))$$(x_{n};f(x_{n}))$  - параллельна оси $OX$$OX$  и при $f(x_n)\ne <!--\; \ne \; -->0$$f(x_{n})\neq 0$  не пересекается с ней в конечной части.
4) Чем больше константа $C>0$$C > 0$  и чем меньше константа $\mathrm{\delta <!--\; \delta \; -->}>0$$\delta >0$  из пункта 3 условий, тем для $x_n\in <!--\; \in \; -->(a,x^{\ast <!--\; \ast \; -->}-<!--\; -\; -->\mathrm{\delta <!--\; \delta \; -->}]\cup <!--\; \cup \; -->[x^{\ast <!--\; \ast \; -->}+\mathrm{\delta <!--\; \delta \; -->},b)$$x_{n}\in (a,\,x^{*}-\delta ]\cup [x^{*}+\delta ,\,b)$  пересечение касательной к графику $f$$f$  и оси $OX$$OX$  ближе к точке $(x^{\ast <!--\; \ast \; -->};0)$$(x^{*};\;0)$  , то есть тем ближе значение $x_{n+1}$$x_{{n+1}}$  к искомой $x^{\ast <!--\; \ast \; -->}\in <!--\; \in \; -->(a,b)$$x^{*}\in (a,\,b)$  .

Итерационный процесс начинается с некоторого начального приближения $x_0\in <!--\; \in \; -->(a,b)$$x_{0}\in (a,\,b)$  , причём между $x_0\in <!--\; \in \; -->(a,b)$$x_{0}\in (a,\,b)$  и искомой точкой $x^{\ast <!--\; \ast \; -->}\in <!--\; \in \; -->(a,b)$$x^{*}\in (a,\,b)$  не должно быть других нулей функции $f$$f$  , то есть «чем ближе $x_0$$x_{0}$  к искомому корню $x^{\ast <!--\; \ast \; -->}$$x^{{*}}$  , тем лучше». Если предположения о нахождении $x^{\ast <!--\; \ast \; -->}$$x^{{*}}$  отсутствуют, методом проб и ошибок можно сузить область возможных значений, применив теорему о промежуточных значениях.

Для предварительно заданных ${\mathrm{\epsilon <!--\; \epsilon \; -->}}_x>0$$\varepsilon _{x}>0$  , ${\mathrm{\epsilon <!--\; \epsilon \; -->}}_f>0$$\varepsilon _{f}>0$  итерационный процесс завершается если $\left\vert {\frac {f(x_{n})}{f'(x_{n})}}\right\vert \approx \vert x_{n+1}-x_{n}\vert <\varepsilon _{x}$  и $\vert f(x_{n+1})\vert <\varepsilon _{f}$  .
В частности, для матрицы дисплея ${\mathrm{\epsilon <!--\; \epsilon \; -->}}_x$$\varepsilon _{x}$  и ${\mathrm{\epsilon <!--\; \epsilon \; -->}}_f$$\varepsilon _{f}$  могут быть рассчитаны, исходя из масштаба отображения графика $f$$f$  , то есть если $x_n$$x_{n}$  и $x_{n+1}$$x_{{n+1}}$  попадают в один вертикальный, а $f(x_n)$$f(x_{n})$  и $f(x_{n+1})$$f(x_{n+1})$  в один горизонтальный ряд.

АлгоритмПравить

1. Задается начальное приближение $x_0$$x_{0}$ .
2. Пока не выполнено условие остановки, в качестве которого можно взять $|x_{n+1}-x_{n}|<\varepsilon$  или $|f(x_{n+1})|<\varepsilon$  (то есть погрешность в нужных пределах), вычисляют новое приближение: $x_{n+1}=x_n-<!--\; -\; -->\frac{f(x_n)}{f^{\prime }(x_n)}$$x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}$ .

ПримерПравить

Иллюстрация применения метода Ньютона к функции $f(x)=\cos <!--\; \; -->x-<!--\; -\; -->x^3$$f(x)=\cos x-x^{3}$  с начальным приближением в точке $x_0=0,5$$x_{0}=0{,}5$ .
График последовательных приближений.	График сходимости.
Согласно способу практического определения скорость сходимости может быть оценена как тангенс угла наклона графика сходимости, то есть в данном случае равна двум.

Рассмотрим задачу о нахождении положительных $x$$x$ , для которых $\cos <!--\; \; -->x=x^3$$\cos x=x^{3}$ . Эта задача может быть представлена как задача нахождения нуля функции $f(x)=\cos <!--\; \; -->x-<!--\; -\; -->x^3$$f(x)=\cos x-x^{3}$ . Имеем выражение для производной $f^{\prime }(x)=-<!--\; -\; -->\sin <!--\; \; -->x-<!--\; -\; -->3x^2$$f'(x)=-\sin x-3x^{2}$ . Так как $\cos <!--\; \; -->x⩽<!--\; ⩽\; -->1$$\cos x\leqslant 1$  для всех $x$$x$  и $x^3>1$$x^{3}>1$  для $x>1$$x>1$ , очевидно, что решение лежит между 0 и 1. Возьмём в качестве начального приближения значение $x_0=0,5$$x_{0}=0{,}5$ , тогда:

${\begin{matrix}x_{1}&=&x_{0}-{\dfrac {f(x_{0})}{f'(x_{0})}}&=&1{,}112\;141\;637\;097,\\x_{2}&=&x_{1}-{\dfrac {f(x_{1})}{f'(x_{1})}}&=&\underline {0{,}}909\;672\;693\;736,\\x_{3}&=&x_{2}-{\dfrac {f(x_{2})}{f'(x_{2})}}&=&\underline {0{,}86}7\;263\;818\;209,\\x_{4}&=&x_{3}-{\dfrac {f(x_{3})}{f'(x_{3})}}&=&\underline {0{,}865\;47}7\;135\;298,\\x_{5}&=&x_{4}-{\dfrac {f(x_{4})}{f'(x_{4})}}&=&\underline {0{,}865\;474\;033\;1}11,\\x_{6}&=&x_{5}-{\dfrac {f(x_{5})}{f'(x_{5})}}&=&\underline {0{,}865\;474\;033\;102}.\end{matrix}}$ 

Подчёркиванием отмечены верные значащие цифры. Видно, что их количество от шага к шагу растёт (приблизительно удваиваясь с каждым шагом): от 1 к 2, от 2 к 5, от 5 к 10, иллюстрируя квадратичную скорость сходимости.

Условия примененияПравить

 

Иллюстрация расхождения метода Ньютона, применённого к функции $f(x)=x^3-<!--\; -\; -->2x+2$$\scriptstyle {f(x)=x^{3}-2x+2}$  с начальным приближением в точке $x_0=0$$\scriptstyle {x_{0}=0}$ .

Рассмотрим ряд примеров, указывающих на недостатки метода.

КонтрпримерыПравить

• Если начальное приближение недостаточно близко к решению, то метод может не сойтись.

Пусть

$f(x)=x^3-<!--\; -\; -->2x+2.$$f(x)=x^{3}-2x+2.$ 

Тогда

$x_{n+1}=x_n-<!--\; -\; -->\frac{x_n^3-<!--\; -\; -->2x_n+2}{3x_n^2-<!--\; -\; -->2}.$$x_{{n+1}}=x_{{n}}-{\frac {x_{n}^{3}-2x_{n}+2}{3x_{n}^{2}-2}}.$ 

Возьмём ноль в качестве начального приближения. Первая итерация даст в качестве приближения единицу. В свою очередь, вторая снова даст ноль. Метод зациклится и решение не будет найдено. В общем случае построение последовательности приближений может быть очень запутанным.

 

График производной функции $f(x)=x+x^2\sin <!--\; \; -->(2/x)$$\scriptstyle {f(x)=x+x^{2}\sin(2/x)}$  при приближении $x$$\scriptstyle {x}$  к нулю справа.

• Если производная не непрерывна в точке корня, то метод может расходиться в любой окрестности корня.

Рассмотрим функцию:

$f(x)={\begin{cases}0,&x=0,\\x+x^{2}\sin \left({\dfrac {2}{x}}\right),&x\neq 0.\end{cases}}$ 

Тогда $f^{\prime }(0)=1$$f'(0)=1$  и $f^{\prime }(x)=1+2x\sin <!--\; \; -->(2/x)-<!--\; -\; -->2\cos <!--\; \; -->(2/x)$$f'(x)=1+2x\sin(2/x)-2\cos(2/x)$  всюду, кроме 0.

В окрестности корня производная меняет знак при приближении $x$$x$  к нулю справа или слева. В то время, как $f(x)⩾<!--\; ⩾\; -->x-<!--\; -\; -->x^2>0$$f(x)\geqslant x-x^{2}>0$  для $0<x<1$ .

Таким образом $f(x)/f^{\prime }(x)$$f(x)/f'(x)$  не ограничено вблизи корня, и метод будет расходиться, хотя функция всюду дифференцируема, её производная не равна нулю в корне, $f$$f$  бесконечно дифференцируема везде, кроме как в корне, а её производная ограничена в окрестности корня.

• Если не существует вторая производная в точке корня, то скорость сходимости метода может быть заметно снижена.

Рассмотрим пример:

$f(x)=x+x^{4/3}.$$f(x)=x+x^{4/3}.$ 

Тогда $f^{\prime }(x)=1+(4/3)x^{1/3}$$f'(x)=1+(4/3)x^{1/3}$  и $f^{″}(x)=(4/9)x^{-<!--\; -\; -->2/3}$$f''(x)=(4/9)x^{-2/3}$  за исключением $x=0$$x=0$ , где она не определена.

На очередном шаге имеем $x_n$$x_{n}$ :

$x_{n+1}=x_n-<!--\; -\; -->\frac{f(x_n)}{f^{\prime }(x_n)}=\frac{(1/3)x_n^{4/3}}{(1+(4/3)x_n^{1/3})}.$$x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}={\frac {(1/3)x_{n}^{4/3}}{(1+(4/3)x_{n}^{1/3})}}.$ 

Скорость сходимости полученной последовательности составляет приблизительно 4/3. Это существенно меньше, нежели 2, необходимое для квадратичной сходимости, поэтому в данном случае можно говорить лишь о линейной сходимости, хотя функция всюду непрерывно дифференцируема, производная в корне не равна нулю, и $f$$f$  бесконечно дифференцируема везде, кроме как в корне.

• Если производная в точке корня равна нулю, то скорость сходимости не будет квадратичной, а сам метод может преждевременно прекратить поиск, и дать неверное для заданной точности приближение.

Пусть

$f(x)=x^2.$$f(x)=x^{2}.$ 

Тогда $f^{\prime }(x)=2x$$f'(x)=2x$  и следовательно $x-<!--\; -\; -->f(x)/f^{\prime }(x)=x/2$$x-f(x)/f'(x)=x/2$ . Таким образом сходимость метода не квадратичная, а линейная, хотя функция всюду бесконечно дифференцируема.

ОграниченияПравить

Пусть задано уравнение $f(x)=0$$f(x)=0$ , где $f(x):<!--\; :\; -->\mathbb{X}\to <!--\; \to \; -->\mathbb{R}$$f(x)\colon \mathbb {X} \to \mathbb {R}$  и надо найти его решение.

Ниже приведена формулировка основной теоремы, которая позволяет дать чёткие условия применимости. Она носит имя советского математика и экономиста Леонида Витальевича Канторовича (1912—1986).

Теорема Канторовича.

Если существуют такие константы $A,B,C$$A,\;B,\;C$ , что:

1. ${\frac {1}{|f'(x)|}}<A$  на $[a,b]$$[a,\;b]$ , то есть $f^{\prime }(x)$$f'(x)$  существует и не равна нулю;
2. $\left|{\frac {f(x)}{f'(x)}}\right|<B$  на $[a,b]$$[a,\;b]$ , то есть $f(x)$$f(x)$  ограничена;
3. $\mathrm{\exists <!--\; \exists \; -->}{f}^{″}(x)$$\exists f''(x)$  на $[a,b]$$[a,\;b]$ , и $|f^{″}(x)|⩽<!--\; ⩽\; -->C⩽<!--\; ⩽\; -->\frac{1}{2AB}$$|f''(x)|\leqslant C\leqslant {\frac {1}{2AB}}$ ;

Причём длина рассматриваемого отрезка $|a-b|<{\frac {1}{AB}}\left(1-{\sqrt {1-2ABC}}\right)$ . Тогда справедливы следующие утверждения:

1. на $[a,b]$$[a,\;b]$  существует корень $x^{\ast <!--\; \ast \; -->}$$x^{*}$  уравнения $f(x)=0:<!--\; :\; -->\mathrm{\exists <!--\; \exists \; -->}{x}^{\ast <!--\; \ast \; -->}\in <!--\; \in \; -->[a,b]:<!--\; :\; -->f(x^{\ast <!--\; \ast \; -->})=0$$f(x)=0\colon \exists x^{*}\in [a,\;b]\colon f(x^{*})=0$ ;
2. если $x_0=\frac{a+b}{2}$$x_{0}={\frac {a+b}{2}}$ , то итерационная последовательность сходится к этому корню: $\left\{x_{n+1}=x_n-<!--\; -\; -->\frac{f(x_n)}{f^{\prime }(x_n)}\right\}\to <!--\; \to \; -->x^{\ast <!--\; \ast \; -->}$$\left\{x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\right\}\to x^{*}$ ;
3. погрешность может быть оценена по формуле $|x^{\ast <!--\; \ast \; -->}-<!--\; -\; -->x_n|⩽<!--\; ⩽\; -->\frac{B}{2^{n-<!--\; -\; -->1}}(2ABC{)}^{2^{n-<!--\; -\; -->1}}$$|x^{*}-x_{n}|\leqslant {\frac {B}{2^{n-1}}}(2ABC)^{2^{n-1}}$ .

Из последнего из утверждений теоремы в частности следует квадратичная сходимость метода:

$|x^{\ast <!--\; \ast \; -->}-<!--\; -\; -->x_n|⩽<!--\; ⩽\; -->\frac{B}{2^{n-<!--\; -\; -->1}}(2ABC{)}^{2^{n-<!--\; -\; -->1}}=\frac{1}{2}\frac{B}{2^{n-<!--\; -\; -->2}}{\left((2ABC{)}^{2^{n-<!--\; -\; -->2}}\right)}^2=\mathrm{\alpha <!--\; \alpha \; -->}|x^{\ast <!--\; \ast \; -->}-<!--\; -\; -->x_{n-<!--\; -\; -->1}{|}^2.$$|x^{*}-x_{n}|\leqslant {\frac {B}{2^{n-1}}}(2ABC)^{2^{n-1}}={\frac {1}{2}}{\frac {B}{2^{n-2}}}\left((2ABC)^{2^{n-2}}\right)^{2}=\alpha |x^{*}-x_{n-1}|^{2}.$ 

Тогда ограничения на исходную функцию $f(x)$$f(x)$  будут выглядеть так:

1. функция должна быть ограничена;
2. функция должна быть гладкой, дважды дифференцируемой;
3. её первая производная $f^{\prime }(x)$$f'(x)$  равномерно отделена от нуля;
4. её вторая производная $f^{″}(x)$$f''(x)$  должна быть равномерно ограничена.

Историческая справкаПравить

Метод был описан Исааком Ньютоном в рукописи «Об анализе уравнениями бесконечных рядов» (лат. «De analysi per aequationes numero terminorum infinitas»), адресованной в 1669 году Барроу, и в работе «Метод флюксий и бесконечные ряды» (лат. «De metodis fluxionum et serierum infinitarum») или «Аналитическая геометрия» (лат. «Geometria analytica») в собраниях трудов Ньютона, которая была написана в 1671 году. В своих работах Ньютон вводит такие понятия, как разложение функции в ряд, бесконечно малые и флюксии (производные в нынешнем понимании). Указанные работы были изданы значительно позднее: первая вышла в свет в 1711 году благодаря Уильяму Джонсону, вторая была издана Джоном Кользоном в 1736 году уже после смерти создателя. Однако описание метода существенно отличалось от его нынешнего изложения: Ньютон применял свой метод исключительно к полиномам. Он вычислял не последовательные приближения $x_n$$x_{n}$ , а последовательность полиномов и в результате получал приближённое решение $x$$x$ .

Впервые метод был опубликован в трактате «Алгебра» Джона Валлиса в 1685 году, по просьбе которого он был кратко описан самим Ньютоном. В 1690 году Джозеф Рафсон опубликовал упрощённое описание в работе «Общий анализ уравнений» (лат. «Analysis aequationum universalis»). Рафсон рассматривал метод Ньютона как чисто алгебраический и ограничил его применение полиномами, однако при этом он описал метод на основе последовательных приближений $x_n$$x_{n}$  вместо более трудной для понимания последовательности полиномов, использованной Ньютоном. Наконец, в 1740 году метод Ньютона был описан Томасом Симпсоном как итеративный метод первого порядка решения нелинейных уравнений с использованием производной в том виде, в котором он излагается здесь. В той же публикации Симпсон обобщил метод на случай системы из двух уравнений и отметил, что метод Ньютона также может быть применён для решения задач оптимизации путём нахождения нуля производной или градиента.

В 1879 году Артур Кэли в работе «Проблема комплексных чисел Ньютона — Фурье» (англ. «The Newton-Fourier imaginary problem») был первым, кто отметил трудности в обобщении метода Ньютона на случай мнимых корней полиномов степени выше второй и комплексных начальных приближений. Эта работа открыла путь к изучению теории фракталов.

Обобщения и модификацииПравить

 

Иллюстрация последовательных приближений метода одной касательной, применённого к функции $f(x)=e^x-<!--\; -\; -->2$$\scriptstyle {f(x)=e^{x}-2}$  с начальным приближением в точке $x_0=1,8$$\scriptstyle {x_{0}=1{,}8}$ .

Метод секущихПравить

Родственный метод секущих является «приближённым» методом Ньютона и позволяет не вычислять производную. Значение производной в итерационной формуле заменяется её оценкой по двум предыдущим точкам итераций:

$f^{\prime }(x_n)\approx <!--\; \approx \; -->\frac{f(x_n)-<!--\; -\; -->f(x_{n-<!--\; -\; -->1})}{x_n-<!--\; -\; -->x_{n-<!--\; -\; -->1}}$$f'(x_{n})\approx {\frac {f(x_{n})-f(x_{n-1})}{x_{n}-x_{n-1}}}$  .

Таким образом, основная формула имеет вид

$x_{n+1}=x_n-<!--\; -\; -->f(x_n)\cdot <!--\; \cdot \; -->\frac{x_n-<!--\; -\; -->x_{n-<!--\; -\; -->1}}{f(x_n)-<!--\; -\; -->f(x_{n-<!--\; -\; -->1})}.$$x_{n+1}=x_{n}-f(x_{n})\cdot {\frac {x_{n}-x_{n-1}}{f(x_{n})-f(x_{n-1})}}.$ 

Этот метод схож с методом Ньютона, но имеет немного меньшую скорость сходимости. Порядок сходимости метода равен золотому сечению — 1,618…

Замечания. 1) Для начала итерационного процесса требуются два различных значения $x_0$$x_{0}$  и $x_1$$x_{1}$  .
2) В отличие от «настоящего метода Ньютона» (метода касательных), требующего хранить только $x_n$$x_{n}$  (и в ходе вычислений — временно $f(x_n)$$f(x_{n})$  и $f^{\prime }(x_n)$$f'(x_{n})$ ), для метода секущих требуется сохранение $x_{n-<!--\; -\; -->1}$$x_{n-1}$  , $x_n$$x_{n}$  , $f(x_{n-<!--\; -\; -->1})$$f(x_{n-1})$  , $f(x_n)$$f(x_{n})$  .
3) Применяется, если вычисление $f^{\prime }(x)$$f'(x)$  затруднено (например, требует большого количества машинных ресурсов: времени и/или памяти).

Метод одной касательнойПравить

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

Формула итераций этого метода имеет вид:

$x_{n+1}=x_n-<!--\; -\; -->\frac{1}{f^{\prime }(x_0)}f(x_n).$$x_{n+1}=x_{n}-{\frac {1}{f'(x_{0})}}f(x_{n}).$ 

Суть метода заключается в том, чтобы вычислять производную лишь один раз, в точке начального приближения $x_0$$x_{0}$ , а затем использовать это значение на каждой последующей итерации:

$\mathrm{\alpha <!--\; \alpha \; -->}(x)={\mathrm{\alpha <!--\; \alpha \; -->}}_0=-<!--\; -\; -->\frac{1}{f^{\prime }(x_0)}.$$\alpha (x)=\alpha _{0}=-{\dfrac {1}{f'(x_{0})}}.$ 

При таком выборе ${\mathrm{\alpha <!--\; \alpha \; -->}}_0$$\alpha _{0}$  в точке $x_0$$x_{0}$  выполнено равенство:

${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x_0)=1+{\mathrm{\alpha <!--\; \alpha \; -->}}_0{f}^{\prime }(x_0)=0,$$\varphi '(x_{0})=1+\alpha _{0}f'(x_{0})=0,$ 

и если отрезок, на котором предполагается наличие корня $x^{\ast <!--\; \ast \; -->}$$x^{*}$  и выбрано начальное приближение $x_0$$x_{0}$ , достаточно мал, а производная ${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x)$$\varphi '(x)$  непрерывна, то значение ${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x^{\ast <!--\; \ast \; -->})$$\varphi '(x^{*})$  будет не сильно отличаться от ${\mathrm{\phi <!--\; \phi \; -->}}^{\prime }(x_0)=0$$\varphi '(x_{0})=0$  и, следовательно, график $y=\mathrm{\phi <!--\; \phi \; -->}(x)$$y=\varphi (x)$  пройдёт почти горизонтально, пересекая прямую $y=x$$y=x$ , что в свою очередь обеспечит быструю сходимость последовательности точек приближений к корню.

Этот метод является частным случаем метода простой итерации. Он имеет линейный порядок сходимости.

Многомерный случайПравить

Обобщим полученный результат на многомерный случай.

Пусть необходимо найти решение системы:

$\left\{{\begin{array}{lcr}f_{1}(x_{1},\;x_{2},\;\ldots ,\;x_{n})&=&0,\\\ldots &&\\f_{m}(x_{1},\;x_{2},\;\ldots ,\;x_{n})&=&0.\end{array}}\right.$ 

Выбирая некоторое начальное значение ${\stackrel{\to <!--\; \to \; -->}{x}}^{[0]}$${\vec {x}}^{[0]}$ , последовательные приближения ${\stackrel{\to <!--\; \to \; -->}{x}}^{[j+1]}$${\vec {x}}^{[j+1]}$  находят путём решения систем уравнений:

$f_i+\underset{k=1}{\overset{n}{\sum <!--\; \sum \; -->}}\frac{\mathrm{\partial <!--\; \partial \; -->}{f}_i}{\mathrm{\partial <!--\; \partial \; -->}{x}_k}(x_k^{[j+1]}-<!--\; -\; -->x_k^{[j]})=0,i=1,2,\dots <!--\; \dots \; -->,m,$$f_{i}+\sum _{k=1}^{n}{\frac {\partial f_{i}}{\partial x_{k}}}(x_{k}^{[j+1]}-x_{k}^{[j]})=0,\qquad i=1,\;2,\;\ldots ,\;m,$ 

где ${\stackrel{\to <!--\; \to \; -->}{x}}^{[j]}=(x_1^{[j]},x_2^{[j]},\dots <!--\; \dots \; -->,x_n^{[j]}),j=0,1,2,\dots <!--\; \dots \; -->$${\vec {x}}^{[j]}=(x_{1}^{[j]},\;x_{2}^{[j]},\;\ldots ,\;x_{n}^{[j]}),\quad j=0,\;1,\;2,\;\ldots$ .

Применительно к задачам оптимизацииПравить

Пусть необходимо найти минимум функции многих переменных $f(\stackrel{\to <!--\; \to \; -->}{x}):<!--\; :\; -->{\mathbb{R}}^n\to <!--\; \to \; -->\mathbb{R}$$f({\vec {x}})\colon \mathbb {R} ^{n}\to \mathbb {R}$ . Эта задача равносильна задаче нахождения нуля градиента $\mathrm{\nabla <!--\; \nabla \; -->}f(\stackrel{\to <!--\; \to \; -->}{x})$$\nabla f({\vec {x}})$ . Применим изложенный выше метод Ньютона:

$\mathrm{\nabla <!--\; \nabla \; -->}f({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})+H({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})({\stackrel{\to <!--\; \to \; -->}{x}}^{[j+1]}-<!--\; -\; -->{\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})=0,j=1,2,\dots <!--\; \dots \; -->,n,$$\nabla f({\vec {x}}^{[j]})+H({\vec {x}}^{[j]})({\vec {x}}^{[j+1]}-{\vec {x}}^{[j]})=0,\quad j=1,\;2,\;\ldots ,\;n,$ 

где $H(\stackrel{\to <!--\; \to \; -->}{x})$$H({\vec {x}})$  — гессиан функции $f(\stackrel{\to <!--\; \to \; -->}{x})$$f({\vec {x}})$ .

В более удобном итеративном виде это выражение выглядит так:

${\stackrel{\to <!--\; \to \; -->}{x}}^{[j+1]}={\stackrel{\to <!--\; \to \; -->}{x}}^{[j]}-<!--\; -\; -->H^{-<!--\; -\; -->1}({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})\mathrm{\nabla <!--\; \nabla \; -->}f({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]}).$${\vec {x}}^{[j+1]}={\vec {x}}^{[j]}-H^{-1}({\vec {x}}^{[j]})\nabla f({\vec {x}}^{[j]}).$ 

Следует отметить, что в случае квадратичной функции метод Ньютона находит экстремум за одну итерацию.

Нахождение матрицы Гессе связано с большими вычислительными затратами, и зачастую не представляется возможным. В таких случаях альтернативой могут служить квазиньютоновские методы, в которых приближение матрицы Гессе строится в процессе накопления информации о кривизне функции.

Метод Ньютона — РафсонаПравить

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

${\stackrel{\to <!--\; \to \; -->}{x}}^{[j+1]}={\stackrel{\to <!--\; \to \; -->}{x}}^{[j]}-<!--\; -\; -->{\mathrm{\lambda <!--\; \lambda \; -->}}_j{H}^{-<!--\; -\; -->1}({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})\mathrm{\nabla <!--\; \nabla \; -->}f({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]}),$${\vec {x}}^{[j+1]}={\vec {x}}^{[j]}-\lambda _{j}H^{-1}({\vec {x}}^{[j]})\nabla f({\vec {x}}^{[j]}),$ 

где ${\mathrm{\lambda <!--\; \lambda \; -->}}_j=\arg <!--\; \; -->\underset{\mathrm{\lambda <!--\; \lambda \; -->}}{min}f({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]}-<!--\; -\; -->\mathrm{\lambda <!--\; \lambda \; -->}{H}^{-<!--\; -\; -->1}({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})\mathrm{\nabla <!--\; \nabla \; -->}f({\stackrel{\to <!--\; \to \; -->}{x}}^{[j]})).$$\lambda _{j}=\arg \min _{\lambda }f({\vec {x}}^{[j]}-\lambda H^{-1}({\vec {x}}^{[j]})\nabla f({\vec {x}}^{[j]})).$  Для оптимизации вычислений применяют следующее улучшение: вместо того, чтобы на каждой итерации заново вычислять гессиан целевой функции, ограничиваются начальным приближением $H(f({\stackrel{\to <!--\; \to \; -->}{x}}^{[0]}))$$H(f({\vec {x}}^{[0]}))$  и обновляют его лишь раз в $m$$m$  шагов, либо не обновляют вовсе.

Применительно к задачам о наименьших квадратахПравить

На практике часто встречаются задачи, в которых требуется произвести настройку свободных параметров объекта или подогнать математическую модель под реальные данные. В этих случаях появляются задачи о наименьших квадратах:

$F(\stackrel{\to <!--\; \to \; -->}{x})=\Vert <!--\; \Vert \; -->\stackrel{\to <!--\; \to \; -->}{f}(\stackrel{\to <!--\; \to \; -->}{x})\Vert <!--\; \Vert \; -->=\underset{i=1}{\overset{m}{\sum <!--\; \sum \; -->}}{f}_i^2(\stackrel{\to <!--\; \to \; -->}{x})=\underset{i=1}{\overset{m}{\sum <!--\; \sum \; -->}}({\mathrm{\phi <!--\; \phi \; -->}}_i(\stackrel{\to <!--\; \to \; -->}{x})-<!--\; -\; -->{\mathcal{F}}_i{)}^2\to <!--\; \to \; -->min.$$F({\vec {x}})=\|{\vec {f}}({\vec {x}})\|=\sum _{{i=1}}^{m}f_{i}^{2}({\vec {x}})=\sum _{{i=1}}^{m}(\varphi _{i}({\vec {x}})-{\mathcal {F}}_{i})^{2}\to \min .$ 

Эти задачи отличаются особым видом градиента и матрицы Гессе:

$\mathrm{\nabla <!--\; \nabla \; -->}F(\stackrel{\to <!--\; \to \; -->}{x})=2J^T(\stackrel{\to <!--\; \to \; -->}{x})\stackrel{\to <!--\; \to \; -->}{f}(\stackrel{\to <!--\; \to \; -->}{x}),$$\nabla F({\vec {x}})=2J^{T}({\vec {x}}){\vec {f}}({\vec {x}}),$ 

$H(\stackrel{\to <!--\; \to \; -->}{x})=2J^T(\stackrel{\to <!--\; \to \; -->}{x})J(\stackrel{\to <!--\; \to \; -->}{x})+2Q(\stackrel{\to <!--\; \to \; -->}{x}),Q(\stackrel{\to <!--\; \to \; -->}{x})=\underset{i=1}{\overset{m}{\sum <!--\; \sum \; -->}}{f}_i(\stackrel{\to <!--\; \to \; -->}{x})H_i(\stackrel{\to <!--\; \to \; -->}{x}),$$H({\vec {x}})=2J^{T}({\vec {x}})J({\vec {x}})+2Q({\vec {x}}),\qquad Q({\vec {x}})=\sum _{{i=1}}^{m}f_{i}({\vec {x}})H_{i}({\vec {x}}),$ 

где $J(\stackrel{\to <!--\; \to \; -->}{x})$$J({\vec {x}})$  — матрица Якоби вектор-функции $\stackrel{\to <!--\; \to \; -->}{f}(\stackrel{\to <!--\; \to \; -->}{x})$${\vec {f}}({\vec {x}})$ , $H_i(\stackrel{\to <!--\; \to \; -->}{x})$$H_{i}({\vec {x}})$  — матрица Гессе для её компоненты $f_i(\stackrel{\to <!--\; \to \; -->}{x})$$f_{i}({\vec {x}})$ .

Тогда очередной шаг $\stackrel{\to <!--\; \to \; -->}{p}$${\vec {p}}$  определяется из системы:

$\left[J^T(\stackrel{\to <!--\; \to \; -->}{x})J(\stackrel{\to <!--\; \to \; -->}{x})+\underset{i=1}{\overset{m}{\sum <!--\; \sum \; -->}}{f}_i(\stackrel{\to <!--\; \to \; -->}{x})H_i(\stackrel{\to <!--\; \to \; -->}{x})\right]\stackrel{\to <!--\; \to \; -->}{p}=-<!--\; -\; -->J^T(\stackrel{\to <!--\; \to \; -->}{x})\stackrel{\to <!--\; \to \; -->}{f}(\stackrel{\to <!--\; \to \; -->}{x}).$$\left[J^{T}({\vec {x}})J({\vec {x}})+\sum _{{i=1}}^{m}f_{i}({\vec {x}})H_{i}({\vec {x}})\right]{\vec {p}}=-J^{T}({\vec {x}}){\vec {f}}({\vec {x}}).$ 

Метод Гаусса — НьютонаПравить

Основная статья: Алгоритм Гаусса — Ньютона

Метод Гаусса — Ньютона строится на предположении о том, что слагаемое $J^T(\stackrel{\to <!--\; \to \; -->}{x})J(\stackrel{\to <!--\; \to \; -->}{x})$$J^{T}({\vec {x}})J({\vec {x}})$  доминирует над $Q(\stackrel{\to <!--\; \to \; -->}{x})$$Q({\vec {x}})$ . Это требование не соблюдается, если минимальные невязки велики, то есть если норма $\Vert <!--\; \Vert \; -->\stackrel{\to <!--\; \to \; -->}{f}(\stackrel{\to <!--\; \to \; -->}{x})\Vert <!--\; \Vert \; -->$$\|{\vec {f}}({\vec {x}})\|$  сравнима с максимальным собственным значением матрицы $J^T(\stackrel{\to <!--\; \to \; -->}{x})J(\stackrel{\to <!--\; \to \; -->}{x})$$J^{T}({\vec {x}})J({\vec {x}})$ . В противном случае можно записать:

$J^T(\stackrel{\to <!--\; \to \; -->}{x})J(\stackrel{\to <!--\; \to \; -->}{x})\stackrel{\to <!--\; \to \; -->}{p}=-<!--\; -\; -->J^T(\stackrel{\to <!--\; \to \; -->}{x})\stackrel{\to <!--\; \to \; -->}{f}(\stackrel{\to <!--\; \to \; -->}{x}).$$J^{T}({\vec {x}})J({\vec {x}}){\vec {p}}=-J^{T}({\vec {x}}){\vec {f}}({\vec {x}}).$ 

Таким образом, когда норма $\Vert <!--\; \Vert \; -->Q(\stackrel{\to <!--\; \to \; -->}{x})\Vert <!--\; \Vert \; -->$$\|Q({\vec {x}})\|$  близка к нулю, а матрица $J(\stackrel{\to <!--\; \to \; -->}{x})$$J({\vec {x}})$  имеет полный столбцевой ранг, шаг $\stackrel{\to <!--\; \to \; -->}{p}$${\vec {p}}$  мало отличается от ньютоновского (с учётом $Q(\stackrel{\to <!--\; \to \; -->}{x})$$Q({\vec {x}})$ ), и метод может достигать квадратичной скорости сходимости, хотя вторые производные и не учитываются. Улучшением метода является алгоритм Левенберга — Марквардта, основанный на эвристических соображениях.

Обобщение на комплексную плоскостьПравить

 

Бассейны Ньютона для полинома пятой степени $p(x)=x^5-<!--\; -\; -->1$$\scriptstyle {p(x)=x^{5}-1}$ . Разными цветами закрашены области притяжения для разных корней. Более тёмные области соответствуют большему числу итераций.

До сих пор в описании метода использовались функции, осуществляющие отображения в пределах множества вещественных значений. Однако метод может быть применён и для нахождения нуля функции комплексной переменной. При этом процедура остаётся неизменной:

$z_{n+1}=z_n-<!--\; -\; -->\frac{f(z_n)}{f^{\prime }(z_n)}.$$z_{n+1}=z_{n}-{\frac {f(z_{n})}{f'(z_{n})}}.$ 

Особый интерес представляет выбор начального приближения $z_0$$z_{0}$ . Ввиду того, что функция может иметь несколько нулей, в различных случаях метод может сходиться к различным значениям, и вполне естественно возникает желание выяснить, какие области обеспечат сходимость к тому или иному корню. Этот вопрос заинтересовал Артура Кэли ещё в 1879 году, однако разрешить его смогли лишь в 70-х годах двадцатого столетия с появлением вычислительной техники. Оказалось, что на пересечениях этих областей (их принято называть областями притяжения) образуются так называемые фракталы — бесконечные самоподобные геометрические фигуры.

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

РеализацияПравить

ScalaПравить

object NewtonMethod {

  val accuracy = 1e-6

  @tailrec
  def method(x0: Double, f: Double => Double, dfdx: Double => Double, e: Double): Double = {
    val x1 = x0 - f(x0) / dfdx(x0)
    if (abs(x1 - x0) < e) x1
    else method(x1, f, dfdx, e)
  }

  def g(C: Double) = (x: Double) => x*x - C

  def dgdx(x: Double) = 2*x

  def sqrt(x: Double) = x match {
    case 0 => 0
    case x if (x < 0) => Double.NaN
    case x if (x > 0) => method(x/2, g(x), dgdx, accuracy) 
  }
}

PythonПравить

from math import sin, cos
from typing import Callable
import unittest


def newton(f: Callable[[float], float], f_prime: Callable[[float], float], x0: float, 
	eps: float=1e-7, kmax: int=1e3) -> float:
	"""
	solves f(x) = 0 by Newton's method with precision eps
	:param f: f
	:param f_prime: f'
	:param x0: starting point
	:param eps: precision wanted
	:return: root of f(x) = 0
	"""
	x, x_prev, i = x0, x0 + 2 * eps, 0
	
	while abs(x - x_prev) >= eps and i < kmax:
		x, x_prev, i = x - f(x) / f_prime(x), x, i + 1

	return x


class TestNewton(unittest.TestCase):
	def test_0(self):
		def f(x: float) -> float:
			return x**2 - 20 * sin(x)


		def f_prime(x: float) -> float:
			return 2 * x - 20 * cos(x)


		x0, x_star = 2, 2.7529466338187049383

		self.assertAlmostEqual(newton(f, f_prime, x0), x_star)


if __name__ == '__main__':
	unittest.main()

PHPПравить

<?php
// PHP 5.4
function newtons_method(
	$a = -1, $b = 1, 
	$f = function($x) {
	
		return pow($x, 4) - 1;
	
	},
	$derivative_f = function($x) {

		return 4 * pow($x, 3);
	
	}, $eps = 1E-3) {

        $xa = $a;
        $xb = $b;

        $iteration = 0;

        while (abs($xb) > $eps) {

            $p1 = $f($xa);
            $q1 = $derivative_f($xa);
            $xa -= $p1 / $q1;
            $xb = $p1;
            ++$iteration;

        }

        return $xa;

}

OctaveПравить

function res = nt()
  eps = 1e-7;
  x0_1 = [-0.5,0.5];
  max_iter = 500;
  xopt = new(@resh, eps, max_iter);   
  xopt
endfunction
function a = new(f, eps, max_iter)
  x=-1;
  p0=1;
  i=0;
 while (abs(p0)>=eps)
    [p1,q1]=f(x);
    x=x-p1/q1;
   p0=p1;
   i=i+1;
 end
 i
 a=x;
endfunction
function[p,q]= resh(x)   % p= -5*x.^5+4*x.^4-12*x.^3+11*x.^2-2*x+1;
   p=-25*x.^4+16*x.^3-36*x.^2+22*x-2;
   q=-100*x.^3+48*x.^2-72*x+22;
endfunction

DelphiПравить

// вычисляемая функция
function fx(x: Double): Double;
begin
  Result := x * x - 17;
end;

// производная функция от f(x)
function dfx(x: Double): Double;
begin
  Result := 2 * x;
end;

function solve(fx, dfx: TFunc<Double, Double>; x0: Double): Double;
const
  eps = 0.000001;
var
  x1: Double;
begin
  x1 := x0 - fx(x0) / dfx(x0); // первое приближение
  while (Abs(x1-x0) > eps) do begin // пока не достигнута точность 0.000001
    x0 := x1;
    x1 := x1 - fx(x1) / dfx(x1); // последующие приближения
  end;
  Result := x1;
end;

// Вызов
solve(fx, dfx,4));

С++Править

#include <stdio.h>
#include <math.h>

#define eps 0.000001
double fx(double x) { return x * x - 17;} // вычисляемая функция
double dfx(double x) { return 2 * x;} // производная функции

typedef double(*function)(double x); // задание типа function

double solve(function fx, function dfx, double x0) {
  double x1  = x0 - fx(x0) / dfx(x0); // первое приближение
  while (fabs(x1 - x0) > eps) { // пока не достигнута точность 0.000001
    x0 = x1;
    x1 = x0 - fx(x0) / dfx(x0); // последующие приближения
  }
  return x1;
}

int main () {
  printf("%f\n", solve(fx, dfx, 4)); // вывод на экран
  return 0;
}

CПравить

typedef double (*function)(double x);

double TangentsMethod(function f, function df, double xn, double eps) {
   double x1  = xn - f(xn)/df(xn);
   double x0 = xn;
   while(abs(x0-x1) > eps) {
      x0 = x1;
      x1 = x1 - f(x1)/df(x1);
   }
   return x1;
}

//Выбор начального приближения
xn = MyFunction(A)*My2Derivative(A) > 0 ? B : A;

double MyFunction(double x) { return (pow(x, 5) - x - 0.2); } //Ваша функция
double MyDerivative(double x) { return (5*pow(x, 4) - 1); } //Первая производная
double My2Derivative(double x) { return (20*pow(x, 3)); } //Вторая производная

//Пример вызова функции
double x = TangentsMethod(MyFunction, MyDerivative, xn, 0.1)

HaskellПравить

import Data.List ( iterate' )

main :: IO ()
main = print $ solve (\ x -> x * x - 17) ( * 2) 4

-- Функция solve универсальна для всех вещественных типов значения которых можно сравнивать.
solve = esolve 0.000001

esolve epsilon func deriv x0 = fst . head $ dropWhile pred pairs
  where
    pred (xn, xn1) = (abs $ xn - xn1) > epsilon -- Функция pred определяет достигнута ли необходимая точность.
    next xn = xn - func xn / deriv xn -- Функция next вычисляет новое приближение.
    iters   = iterate' next x0        -- Бесконечный список итераций.
    pairs   = zip iters (tail iters)  -- Бесконечный список пар итераций вида: [(x0, x1), (x1, x2) ..].

ЛитератураПравить

• Акулич И. Л. Математическое программирование в примерах и задачах : Учеб. пособие для студентов эконом. спец. вузов. — М. : Высшая школа, 1986. — 319 с. : ил. — ББК 22.1 А44. — УДК 517.8.
• Амосов А. А., Дубинский Ю. А., Копченова Н. П. Вычислительные методы для инженеров : Учеб. пособие. — М. : Высшая школа, 1994. — 544 с. : ил. — ББК 32.97 А62. — УДК 683.1. — ISBN 5-06-000625-5.
• Бахвалов Н. С., Жидков Н. П., Кобельков Г. Г. Численные методы. — 8-е изд. — М. : Лаборатория Базовых Знаний, 2000.
• Вавилов С. И. Исаак Ньютон. — М. : Изд. АН СССР, 1945.
• Волков Е. А. Численные методы. — М. : Физматлит, 2003.
• Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. Пер. с англ. — М. : Мир, 1985.
• Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. — М. : Наука, 1970. — С. 575—576.
• Коршунов Ю. М., Коршунов Ю. М. Математические основы кибернетики. — Энергоатомиздат, 1972.
• Максимов Ю. А.,Филлиповская Е. А. Алгоритмы решения задач нелинейного программирования. — М. : МИФИ, 1982.
• Морозов А. Д. Введение в теорию фракталов. — МИФИ, 2002.

См. такжеПравить

• Метод простой итерации
• Двоичный поиск
• Метод дихотомии
• Метод секущих
• Метод Мюллера
• Метод хорд
• Метод бисекции
• Фрактал
• Численное решение уравнений
• Целочисленный квадратный корень
• Быстрый обратный квадратный корень

СсылкиПравить

• «Бассейны Ньютона» на сайте fractalworld.xaoc.ru
• «Исаак Ньютон» на сайте www.scottish-wetlands.org
• «Математические работы Канторовича» на сайте Института математики СО РАН
• Hazewinkel, Michiel, ed. (2001), Newton method, Encyclopedia of Mathematics, Springer, ISBN 978-1-55608-010-4 
• Weisstein, Eric W. Newton's Method (англ.) на сайте Wolfram MathWorld.
• Newton’s method, Citizendium.
• Mathews, J., The Accelerated and Modified Newton Methods, Course notes.
• Wu, X., Roots of Equations, Course notes.