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

Материал из Википедии — свободной энциклопедии
Перейти к: навигация, поиск

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

Описание метода[править | править вики-текст]

Обоснование[править | править вики-текст]

Чтобы численно решить уравнение f(x)=0\! методом простой итерации, его необходимо привести к следующей форме: x=\varphi(x)\!, где \varphi\! — сжимающее отображение.

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

\varphi'(x^*)=1+\alpha'(x^*)f(x^*)+\alpha(x^*) f'(x^*)=0.\!

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

\alpha(x)=-\frac{1}{f'(x)}.\!

С учётом этого функция \varphi(x)\! определяется:

\varphi(x)=x-\frac{f(x)}{f'(x)}.\!

При некоторых условиях эта функция в окрестности корня осуществляет сжимающее отображение[1], и алгоритм нахождения численного решения уравнения f(x)=0\! сводится к итерационной процедуре вычисления:

x_{n+1}=x_{n}-\frac{f(x_n)}{f'(x_n)}.\!

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

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

Геометрическая интерпретация[править | править вики-текст]

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

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

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}, \, \!

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

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

x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}. \, \!

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

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

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

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

Для предварительно заданных \varepsilon_{x} >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}\! .
В частности, для матрицы дисплея \varepsilon_{x} \! и \varepsilon_{f} \! могут быть рассчитаны, исходя из масштаба отображения графика f , т.е. если x_n и x_{n+1} - попадают в один вертикальный, а f(x_n) и f(x_{n+1}) - в один горизонтальный ряд.

Алгоритм[править | править вики-текст]

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

Пример[править | править вики-текст]

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

Рассмотрим задачу о нахождении положительных x, для которых \cos x=x^3. Эта задача может быть представлена как задача нахождения нуля функции f(x)=\cos x-x^3. Имеем выражение для производной f'(x)=-\sin x-3x^2. Так как \cos x\leqslant 1 для всех x и x^3>1 для x>1, очевидно, что решение лежит между 0 и 1. Возьмём в качестве начального приближения значение 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, иллюстрируя квадратичную скорость сходимости.


Условия применения[править | править вики-текст]

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

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

Контрпримеры[править | править вики-текст]

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

Пусть

f(x)=x^3-2x+2.\!

Тогда

x_{n+1}=x_{n}-\frac{x_n^3-2x_n+2}{3x_n^2-2}.

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

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

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

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

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

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

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

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

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

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

На очередном шаге имеем x_n\!:

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(x)=x^2.\!

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

Ограничения[править | править вики-текст]

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

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

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

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

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

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

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

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

|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)\! будут выглядеть так:

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

Историческая справка[править | править вики-текст]

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

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

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

Обобщения и модификации[править | править вики-текст]

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

Метод секущих[править | править вики-текст]

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

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 \frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})} . \, \!

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

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

Метод одной касательной[править | править вики-текст]

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

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

x_{n+1}=x_n-\frac{1}{f'(x_0)}f(x_n).\!

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

\alpha(x)=\alpha_0=-\dfrac{1}{f'(x_0)}.\!

При таком выборе \alpha_0\! в точке x_0\! выполнено равенство:

\varphi'(x_0)=1+\alpha_0f'(x_0)=0,\!

и если отрезок, на котором предполагается наличие корня x^*\! и выбрано начальное приближение x_0\!, достаточно мал, а производная \varphi'(x)\! непрерывна, то значение \varphi'(x^*)\! будет не сильно отличаться от \varphi'(x_0)=0\!\! и, следовательно, график  y=\varphi(x)\! пройдёт почти горизонтально, пересекая прямую y=x\!, что в свою очередь обеспечит быструю сходимость последовательности точек приближений к корню.

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

Многомерный случай[править | править вики-текст]

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

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

\left\{\begin{array}{lcr}
f_1(x_1,\;x_2,\;\ldots,\;x_n) & = & 0, \\
\ldots & & \\
f_n(x_1,\;x_2,\;\ldots,\;x_n) & = & 0.
\end{array}\right.

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

f_i+\sum_{k=1}^n\frac{\partial f_i}{\partial x_k}(x^{[j+1]}_k-x_k^{[j]})=0,\qquad i=1,\;2,\;\ldots,\;n,\!

где \vec{x}^{[j]}=(x_1^{[j]},\;x_2^{[j]},\;\ldots,\;x_n^{[j]}),\quad j=0,\;1,\;2,\;\ldots\!.

Применительно к задачам оптимизации[править | править вики-текст]

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

\nabla f(\vec{x}^{[j]})+H(\vec{x}^{[j]})(\vec{x}^{[j+1]}-\vec{x}^{[j]})=0,\quad j=1,\;2,\;\ldots,\;n,\!

где H(\vec{x})\! — гессиан функции f(\vec{x})\!.

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

\vec{x}^{[j+1]}=\vec{x}^{[j]}-H^{-1}(\vec{x}^{[j]})\nabla f(\vec{x}^{[j]}).\!

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

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

Метод Ньютона — Рафсона[править | править вики-текст]

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

\vec{x}^{[j+1]}=\vec{x}^{[j]}-\lambda_j H^{-1}(\vec{x}^{[j]})\nabla f(\vec{x}^{[j]}),\!

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

Применительно к задачам о наименьших квадратах[править | править вики-текст]

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

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.

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

\nabla F(\vec{x})=2J^T(\vec{x})\vec{f}(\vec{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(\vec{x}) — матрица Якоби вектор-функции \vec{f}(\vec{x}), H_i(\vec{x}) — матрица Гессе для её компоненты f_i(\vec{x}).

Тогда очередное направление \vec{p} определяется из системы:

\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(\vec{x})J(\vec{x}) доминирует над Q(\vec{x}). Это требование не соблюдается, если минимальные невязки велики, то есть если норма \|\vec{f}(\vec{x})\| сравнима с максимальным собственным значением матрицы J^T(\vec{x})J(\vec{x}). В противном случае можно записать:

J^T(\vec{x})J(\vec{x})\vec{p}=-J^T(\vec{x})\vec{f}(\vec{x}).

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

Обобщение на комплексную плоскость[править | править вики-текст]

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

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

z_{n+1}=z_n-\frac{f(z_n)}{f'(z_n)}.\!

Особый интерес представляет выбор начального приближения 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[править | править вики-текст]

def newtons_method(x0, f, f1, e):
    #f1 - производная
    while True:
        x1 = x0 - (f(x0) / f1(x0))
        if abs(x1 - x0) < e:
            return x1
        x0 = x1

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

С[править | править вики-текст]

#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);

double solve(function fx, function dfx, double x0) {
  double x1  = x0 - fx(x0)/dfx(x0);
  while (fabs(x1-x0)>eps) {
    x0 = x1;
    x1 = x1 - fx(x1)/dfx(x1);
  }
  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)

Литература[править | править вики-текст]

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

Примечания[править | править вики-текст]

  1. Доказательство: Пусть дана функция вещественного переменного дважды непрерывно дифференцируемая в своей области определения, производная которой нигде не обращается в нуль:
    \scriptstyle{f(x)\colon\mathbb{X}\to\R,\;f(x)\in\mathrm{C}^2(\mathbb{X});\quad\forall x\in\mathbb{X}\;f'(x)\neq 0.}
    И необходимо доказать, что функция \scriptstyle{\varphi(x)=x-\frac{f(x)}{f'(x)}} осуществляет сжимающее отображение вблизи корня уравнения \scriptstyle{f(x)=0}. В силу непрерывной дифференцируемости функции \scriptstyle{f(x)} и неравенства нулю её первой производной \scriptstyle{\varphi(x)} непрерывна. Производная \scriptstyle{\varphi'(x)} равна:
    \scriptstyle{\varphi'(x)=\frac{f(x)f''(x)}{\left(f'(x)\right)^2}.}
    В условиях, наложенных на \scriptstyle{f(x)}, она также непрерывна. Пусть \scriptstyle{\tilde{x}} — искомый корень уравнения: \scriptstyle{f(\tilde{x})=0}, следовательно в его окрестности \scriptstyle{\varphi'(x)\approx 0}:
    \scriptstyle{\forall\varepsilon\colon 0<\varepsilon<1,\;\exist\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})\;\exist\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|.}
    В силу того, что \scriptstyle{\varphi(\tilde{x})=\tilde{x}} в этой же дельта окрестности выполняется:
    \scriptstyle{\forall x\in U_{\delta}(\tilde{x})\colon\;|\varphi(x)-\tilde{x}|<\varepsilon|x-\tilde{x}|.}
    Таким образом полученная функция \scriptstyle{\varphi(x)} в окрестности корня \scriptstyle{U_\delta(\tilde{x})} осуществляет сжимающее отображение.

См. также[править | править вики-текст]

Ссылки[править | править вики-текст]