Участник:Gukov/Песочница
Материал из MachineLearning.
(→Числовой пример) |
(→Случай почти равных корней) |
||
(23 промежуточные версии не показаны) | |||
Строка 1: | Строка 1: | ||
== Введение == | == Введение == | ||
- | === | + | === Метод Ньютона-Рафсона === |
+ | {{main|Метод касательных (Ньютона-Рафсона)}} | ||
- | + | Пусть имеется некоторая функция <tex>F(x)</tex> и необходимо найти такие значения аргумента <tex>x</tex>, для которых | |
- | {{ eqno | 1}} | + | {{eqno|1}} |
- | :<tex> | + | :<tex>F(x)=0</tex> |
- | + | Перепишем {{eqref|1}} в виде | |
- | <tex> | + | {{eqno|1.1}} |
- | f(x) | + | :<tex>x = f(x)</tex> |
- | </tex> | + | |
- | {{ eqno | 2}} | + | и запишем <tex>n+1</tex>-ое приближения корня {{eqref|1.1}}, при этом делая поправку <tex>\alpha</tex> к очередному значению <tex>x_n</tex> |
- | :<tex> | + | {{eqno|2}} |
+ | :<tex>x_{n+1} = x_n + \alpha \Delta x_n,</tex> | ||
+ | где <tex>\Delta x_n = f(x_n) - x_n,</tex> и положим <tex>\alpha = \frac{1}{1-f'(x_n)}</tex>. | ||
- | + | Тогда {{eqref|2}} перепишется в виде | |
- | + | {{eqno|3}} | |
+ | :<tex>x_{n+1} = \frac{f(x_n) - x_n f'(x_n)}{1 - f'(x_n)}</tex> | ||
- | :<tex> | + | Нетрудно видеть, что {{eqref|3}} эквивалентно простому [[Метод последовательных приближений|методу последовательных приближений]] |
+ | :<tex>x_{n+1} = g(x_n),</tex> | ||
+ | где <tex>g(x) = \frac{f(x) - x f'(x)}{1 - f'(x)}</tex>. | ||
- | + | Вспомним также, что если <tex>|g'(x)| < 1</tex>, то метод сходится. Имеем | |
- | + | :<tex>g'(x) = \frac{f''(x)[f(x) - x]}{[1 - f'(x)]^2}.</tex> | |
- | + | Но так как справедливо соотношение {{eqref|1.1}}, то для <tex>x,</tex> достаточно близких к решению {{eqref|1.1}}, выражение в скобках в числителе дроби становится малым. Поэтому итерационный метод, описываемый формулой {{eqref|3}} сходится, если | |
+ | * 1. Начальное приближение <tex>x_0</tex> выбрано достаточно близко к решению <tex>x = f(x)</tex>. | ||
+ | * 2. Производная <tex>f''(x)</tex> не становится слишком большой. | ||
+ | * 3. Производная <tex>f'(x)</tex> не слишком близка к 1. | ||
- | + | Это и есть метод Ньютона-Рафсона. Обычно его записывают в виде | |
+ | {{eqno|4}} | ||
+ | :<tex>x_{n+1} = x_{n} - \frac{F(x_n)}{F'(x_n)}</tex>, | ||
+ | |||
+ | где <tex>F(x) = f(x) - x = 0</tex> | ||
+ | |||
+ | Таким образом, мы вернулись к уравнению в форме {{eqref|1}}, и условия сходимости принимают следующий вид | ||
+ | * 1. Начальное приближение <tex>x_0</tex> выбрано достаточно близко к корню уравнения <tex>F(x) = 0</tex>. | ||
+ | * 2. Производная <tex>F''(x)</tex> не становится очень большой. | ||
+ | * 3. Производная <tex>F'(x)</tex> не слишком близка к 0. | ||
+ | |||
+ | === Случай почти равных корней === | ||
+ | [[Изображение:macon.png|thumb|200px|Рисунок 1]] | ||
+ | Условие 3 сходимости метода Ньютона-Рафсона означает, что никакие два корня не находятся слишком близко один к другому. Соответствующая ситуация представлена на рисунке 1 (масштаб сильно увеличен). Заметим, что производная <tex>f'(x)</tex> близка | ||
+ | к 1 при x, равном обоим значениям корней, <tex>a_1</tex> и <tex>a_2</tex>. Более того, на основании теоремы Лагранжа, можно утверждать, что <tex>f'(x) = 1</tex> где-то между <tex>a_1</tex> и <tex>a_2</tex>. | ||
+ | |||
+ | Рассмотрим, что случится, если принять <tex>x_0</tex> в качестве исходного значения для корня <tex>a_1</tex>. Касательная, | ||
+ | проведенная через точку <tex>C</tex>, пересечет прямую <tex>y=x</tex> в точке <tex>A</tex>, и следующее приближение будет равно | ||
+ | <tex>x_1</tex>. Касательная, проведенная через точку <tex>B</tex>, пересекает прямую в точке <tex>D</tex>, и в качестве следующего приближеня получается снова <tex>x_0</tex>. Итерационный процесс, таким образом, осциллирует между <tex>x_0</tex> и <tex>x_1</tex> | ||
+ | до бесконечности, не сходясь ни к одному значению корня. Иначе говоря, не удается <i>отделить</i> эти два корня, потому что они расположены слишком близко. | ||
+ | |||
+ | Поэтому необходимо начальное приближение, достаточно близкое к искомому значению корня. Трудности возникают потому, что | ||
+ | вычисление знаменателя в формуле {{eqref|3}} включает в себя вычитание двух почти равных чисел, что приводит к понижению точности. | ||
== Изложение метода == | == Изложение метода == | ||
- | === | + | === Поиск начального приближения === |
- | + | [[Изображение:macon2.PNG|thumb|200px|Рисунок 2]] | |
- | <tex> | + | Сначала находится значение x, при котором <tex>f'(x)=1</tex>, то есть решается уравнение |
- | + | :<tex>x=x+f'(x)-1</tex> | |
- | + | Пусть решением этого уравнения будет некоторое <tex>x = x_0</tex>. Эта точка расположена между двумя корнями, <tex>a_1</tex> и <tex>a_2</tex>. Чтобы получить начальное приближение для решения уравнения, предположим, что <tex>x_0</tex> лежит посредине между | |
- | что | + | <tex>a_1</tex> и <tex>a_2</tex> (рисунок 2). Другими словами, мы предполагаем, что <tex>x_0 + d</tex> и <tex>x_0 - d</tex> являются корнями уравнения {{eqref|1.1}}. Разлагая <tex>f(x)</tex> в ряд Тейлора в окрестности точки <tex>x_0</tex> и принимая во внимание, что <tex>f'(x)=1</tex>, получаем |
+ | :<tex>f(x) = f(x_0) + (x - x_0) + \frac{1}{2} f''(x_0) (x - x_0)^2 + \ldots</tex> | ||
- | + | Ограничим ряд тремя членами. Подставляя <tex>x + d</tex> вместо <tex>x</tex>, имеем | |
- | :<tex> | + | :<tex>f(x_0+d) = f(x_0) + d + \frac{1}{2} f''(x_0)(d^2)</tex> |
- | + | Но по условию | |
- | + | :<tex>f(x_0+d) = x_0+d,</tex> | |
- | + | ||
- | :<tex> | + | |
- | + | поэтому, решая эти уравнения относительно <tex>d</tex>, получаем | |
+ | :<tex>d=\sqrt{\frac{2(x_0 - f(x_0))}{f''(x_0)}}</tex> | ||
- | :<tex> | + | Таким образом, процесс решения сводится к следующему. Если дано уравнение с почти равными корнями, то, определив приблизительное |
+ | местонахождение этих корней, необходиом решить уравнение | ||
+ | :<tex>x = x + f'(x) - 1</tex> | ||
- | + | и определить значение <tex>x_0</tex>. Для решения этого уравнения можно применить, например, метод Ньютона-Рафсона. | |
+ | Найдя значение <tex>x_0</tex>, можно определить значение <tex>d</tex>. И, наконец, значения <tex>x_0-d</tex> и <tex>x_0+d</tex> | ||
+ | используются в качестве начальных приближений для определения соответственно <tex>a_1</tex> и <tex>a_2</tex>. | ||
- | + | === Метод Рунге === | |
- | + | Пусть существует асимптотическое разложение вида: | |
- | + | <tex>D_N(f) = I_N(f) - I_0 = \alpha _2 h^2 + \alpha _4 h^4 + \alpha _6 h^6 + \ldots</tex>, | |
- | + | где <tex>f(x)</tex> - достаточно гладкая функция, а <tex>hN = b - a</tex>. | |
- | + | ||
- | + | Проведем расчеты на двух равномерных сетках с шагами <tex>h_1</tex> и <tex>h_2</tex> соответственно и найдем выражения | |
+ | <tex>I^{h_1}[f] = I_{N_1}[f]</tex> и <tex>I^{h_2}[f] = I_{N_2}[f]</tex>, <tex>h_1 N_1 = h_2 N_2 = b - a</tex>. Потребуем, | ||
+ | чтобы погрешность для их линейной комбинации: | ||
- | + | :<tex>{\tilde D}^{h}(f) = \sigma D^{h_1}(f) + (1 - \sigma)D^{h_2}(f)</tex> | |
- | < | + | была величиной более высокого порядка по сравнению с <tex>D^{h_1}</tex> и <tex>D^{h_2}</tex>. Если для <tex>D^h = D_N</tex> имеет место формула вида |
- | <tex> | + | |
- | < | + | :<tex> D^{h} = I^{h}[f] - I[f] = \alpha _p h^p + \alpha _q h^q + \ldots,</tex> <tex>q > p</tex>, |
- | + | ||
- | + | то для | |
- | + | ||
- | + | :<tex>{\tilde D}^{h} = (\sigma I^{h_1}[f] + (1 - \sigma)I^{h_2}[f]) - I[f]</tex> | |
- | + | ||
- | + | получим | |
- | + | ||
- | + | :<tex>{\tilde D}^{h}(f) = \alpha _p(\sigma h_1^{p} + (1 - \sigma)h_2^{p}) + \alpha _q (\sigma h_1^{q} + (1 - \sigma)h_2^{q}) + \ldots</tex> | |
- | + | ||
- | + | Выберем параметр \sigma из условия <tex>\sigma h_1^{p} + (1 - \sigma)h_2^{p}</tex>: | |
- | + | ||
- | + | :<tex>\sigma = \frac{h_2^p}{h_2^p - h_1^p}</tex>. | |
- | + | ||
- | + | Имеем | |
- | + | ||
- | + | :<tex>{\tilde D}^h(f) = \alpha _q(\sigma h_1^q + (1 - \sigma)h_2^q) + \ldots = O(h^q)</tex>, <tex>h = max(h_1, h_2)</tex>, | |
- | + | ||
- | + | причем <tex>\sigma h_1^q + (1 - \sigma)h^q_2 < 0</tex>. Так, если <tex>p = 2, q = 4</tex>, то | |
- | </ | + | <tex>{\tilde D}^{h}(f) = - \alpha _4 h_1^2 h_2^2 + \ldots = O(h^4)</tex>. Таким образом, проведя вычисления на двух сетках |
- | <tex>\ | + | с шагами <tex>h_1</tex> и <tex>h_2 \ne h_1</tex>, мы повысили порядок точности на <tex>2</tex> (на <tex>(q - p)</tex>) для |
- | = | + | <tex>{\tilde I} = \sigma I^{h_1} + (1 - \sigma)I^{h_2}</tex>. |
+ | |||
+ | === Процесс Эйткена === | ||
+ | |||
+ | Метод расчета на нескольких сетках применяется для повышения порядка точности даже в том случае, когда неизвестен порядок главного члена погрешности. Предположим, чт для погрешности имеет место представление | ||
+ | |||
+ | :<tex>D^{h}(f) = \alpha _p h^{p} + \alpha _q h^{q} + \ldots</tex>, <tex>q > p</tex>, | ||
+ | |||
+ | так что | ||
+ | |||
+ | :<tex>I^{h}[f] = I[f] + \alpha _p h^p + \alpha _q h^{q} + \ldots</tex> | ||
+ | |||
+ | Проведем вычисления на трех сетках: <tex>h_1 = h</tex>, <tex>h_2 = \rho h</tex>, <tex>h_3 = \rho^{2} h</tex> (<tex>0 < \rho < 1</tex>). Определим <tex>p</tex>. При этом пренебрегаем членом <tex>O(h^q)</tex>. Образуем отношение | ||
+ | |||
+ | :<tex>A = \frac{I^{h_1}[f] - I^{h_2}[f]}{I^{h_2}[f] - I^{h_3}[f]} \approx \frac{h^p_1 - h^p_2}{h^p_2 - h^p_3} = \frac{1 - \rho^p}{\rho^p(1-\rho^p)} = (\frac{1}{\rho})^p</tex> | ||
+ | |||
+ | и найдем | ||
+ | |||
+ | :<tex>p \approx \frac{\ln A}{\ln(\frac{1}{\rho})}</tex>. | ||
+ | |||
+ | Зная приближенное значение <tex>p</tex>, можно методом методом Рунге повысить порядок точности. Для этого образуем комбинацию | ||
+ | <tex>{\tilde I}^h = \sigma I^{h_1} + (1 - \sigma)I^{h_2}</tex> и выберем <tex>\sigma</tex> так, чтобы <tex> \sigma h_1^p+(1 - \sigma)h_2^p = (\sigma + (1 - \sigma)\rho^p)h^p = 0 </tex>, то есть <tex>\sigma = \frac{\rho^p}{\rho^p - 1} = \frac{1}{1 - A}</tex>. | ||
+ | Тогда для погрешности <tex>{\tilde D}^h = {\tilde I}^h - I</tex> получаем | ||
+ | |||
+ | :<tex> {\tilde D}^h(f) = O(h^q) </tex>. | ||
== Числовой пример == | == Числовой пример == | ||
Строка 95: | Строка 154: | ||
:<tex>\int^{5}_{1} \ln x\, dx = x \ln x - x |_1^{5} = 5\ln 5 - 4</tex> | :<tex>\int^{5}_{1} \ln x\, dx = x \ln x - x |_1^{5} = 5\ln 5 - 4</tex> | ||
+ | [[Изображение:Extrapol.png|thumb|300px]] | ||
В нижеследующей таблице представлены результаты работы программы: | В нижеследующей таблице представлены результаты работы программы: | ||
- | + | {|class="standard" | |
- | + | !r | |
- | + | !Исходная формула | |
- | + | !Экстраполированная формула | |
- | + | !Точное значение | |
+ | !Погрешность вычислений | ||
+ | !Погрешность формулы | ||
|- | |- | ||
- | + | !2 | |
- | | | + | |3.98277278 |
- | | | + | |4.04665506 |
- | | | + | |4.04718956 |
+ | |0.0005345 | ||
+ | |0.00275556 | ||
|- | |- | ||
- | |4 | + | !4 |
- | | | + | |4.03068449 |
- | | | + | |4.04714980 |
- | | | + | |4.04718956 |
+ | |0.00003976 | ||
+ | |0.00017222 | ||
|- | |- | ||
- | + | !8 | |
- | | | + | |4.04303347 |
- | | | + | |4.04718692 |
- | |4. | + | |4.04718956 |
+ | |0.00000264 | ||
+ | |0.00001076 | ||
|- | |- | ||
- | + | !16 | |
- | | | + | |4.04614856 |
- | | | + | |4.04718939 |
- | |4. | + | |4.04718956 |
+ | |0.00000017 | ||
+ | |0.00000067 | ||
|- | |- | ||
- | + | !32 | |
- | | | + | |4.04692918 |
- | | | + | |4.04718955 |
- | |4. | + | |4.04718956 |
+ | |0.00000001 | ||
+ | |0.00000004 | ||
|- | |- | ||
- | + | !64 | |
- | | | + | |4.04712446 |
- | |4. | + | |4.04718956 |
- | |4. | + | |4.04718956 |
+ | |0 | ||
+ | |0 | ||
|- | |- | ||
- | + | !20384 | |
- | + | |4.04718956 | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | |4. | + | |
| | | | ||
- | |||
- | |||
- | |||
- | |||
| | | | ||
- | |||
- | |||
- | |||
| | | | ||
| | | | ||
- | |} | + | |} |
+ | |||
+ | Здесь <tex>r</tex> - коэффициент измельчения шага <tex>h</tex>. Исходная величина шага <tex>h = 2</tex>. | ||
+ | |||
+ | На иллюстрации черная сплошная линия - исходная формула, красная пунктирная - экстраполированная. | ||
+ | |||
+ | Как мы видим, разница между экстраполированными и неэкстраполированными результатами значительна. Уже при величине шага в <tex>h = \frac{1}{64}</tex> мы можем найти значение интеграла с точностью <tex>10^{-8}</tex>, тогда как в исходной формуле нам для достижения такой точности пришлось бы задать величину шага <tex>h = \frac{1}{10192}</tex>. | ||
== Рекомендации программисту == | == Рекомендации программисту == | ||
+ | |||
+ | Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии. | ||
+ | Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца. | ||
== Заключение == | == Заключение == | ||
+ | |||
+ | Мы получили более точные результаты при меньшем количестве вычислений функции, чем в исходном методе. Избавившись от степени, составляющей ошибку, мы пришли к результату, который в противном случае был бы недостижим без значительного уменьшения размера шага. Таким образом, мы преобразовали весьма посредственный алгоритм вычисления интегралов в довольно эффективный. | ||
== Список литературы == | == Список литературы == |
Текущая версия
Содержание |
Введение
Метод Ньютона-Рафсона
Пусть имеется некоторая функция и необходимо найти такие значения аргумента , для которых
Перепишем (1) в виде
и запишем -ое приближения корня (1.1), при этом делая поправку к очередному значению
где и положим .
Тогда (2) перепишется в виде
Нетрудно видеть, что (3) эквивалентно простому методу последовательных приближений
где .
Вспомним также, что если , то метод сходится. Имеем
Но так как справедливо соотношение (1.1), то для достаточно близких к решению (1.1), выражение в скобках в числителе дроби становится малым. Поэтому итерационный метод, описываемый формулой (3) сходится, если
- 1. Начальное приближение выбрано достаточно близко к решению .
- 2. Производная не становится слишком большой.
- 3. Производная не слишком близка к 1.
Это и есть метод Ньютона-Рафсона. Обычно его записывают в виде
- ,
где
Таким образом, мы вернулись к уравнению в форме (1), и условия сходимости принимают следующий вид
- 1. Начальное приближение выбрано достаточно близко к корню уравнения .
- 2. Производная не становится очень большой.
- 3. Производная не слишком близка к 0.
Случай почти равных корней
Условие 3 сходимости метода Ньютона-Рафсона означает, что никакие два корня не находятся слишком близко один к другому. Соответствующая ситуация представлена на рисунке 1 (масштаб сильно увеличен). Заметим, что производная близка к 1 при x, равном обоим значениям корней, и . Более того, на основании теоремы Лагранжа, можно утверждать, что где-то между и .
Рассмотрим, что случится, если принять в качестве исходного значения для корня . Касательная, проведенная через точку , пересечет прямую в точке , и следующее приближение будет равно . Касательная, проведенная через точку , пересекает прямую в точке , и в качестве следующего приближеня получается снова . Итерационный процесс, таким образом, осциллирует между и до бесконечности, не сходясь ни к одному значению корня. Иначе говоря, не удается отделить эти два корня, потому что они расположены слишком близко.
Поэтому необходимо начальное приближение, достаточно близкое к искомому значению корня. Трудности возникают потому, что вычисление знаменателя в формуле (3) включает в себя вычитание двух почти равных чисел, что приводит к понижению точности.
Изложение метода
Поиск начального приближения
Сначала находится значение x, при котором , то есть решается уравнение
Пусть решением этого уравнения будет некоторое . Эта точка расположена между двумя корнями, и . Чтобы получить начальное приближение для решения уравнения, предположим, что лежит посредине между и (рисунок 2). Другими словами, мы предполагаем, что и являются корнями уравнения (1.1). Разлагая в ряд Тейлора в окрестности точки и принимая во внимание, что , получаем
Ограничим ряд тремя членами. Подставляя вместо , имеем
Но по условию
поэтому, решая эти уравнения относительно , получаем
Таким образом, процесс решения сводится к следующему. Если дано уравнение с почти равными корнями, то, определив приблизительное местонахождение этих корней, необходиом решить уравнение
и определить значение . Для решения этого уравнения можно применить, например, метод Ньютона-Рафсона. Найдя значение , можно определить значение . И, наконец, значения и используются в качестве начальных приближений для определения соответственно и .
Метод Рунге
Пусть существует асимптотическое разложение вида:
,
где - достаточно гладкая функция, а .
Проведем расчеты на двух равномерных сетках с шагами и соответственно и найдем выражения и , . Потребуем, чтобы погрешность для их линейной комбинации:
была величиной более высокого порядка по сравнению с и . Если для имеет место формула вида
- ,
то для
получим
Выберем параметр \sigma из условия :
- .
Имеем
- , ,
причем . Так, если , то . Таким образом, проведя вычисления на двух сетках с шагами и , мы повысили порядок точности на (на ) для .
Процесс Эйткена
Метод расчета на нескольких сетках применяется для повышения порядка точности даже в том случае, когда неизвестен порядок главного члена погрешности. Предположим, чт для погрешности имеет место представление
- , ,
так что
Проведем вычисления на трех сетках: , , (). Определим . При этом пренебрегаем членом . Образуем отношение
и найдем
- .
Зная приближенное значение , можно методом методом Рунге повысить порядок точности. Для этого образуем комбинацию и выберем так, чтобы , то есть . Тогда для погрешности получаем
- .
Числовой пример
Найдем с помощью квадратурной формулы трапеций приближенное значение интеграла, применив экстраполяцию Ричардсона (данный метод называется методом Ромберга):
В нижеследующей таблице представлены результаты работы программы:
r | Исходная формула | Экстраполированная формула | Точное значение | Погрешность вычислений | Погрешность формулы |
---|---|---|---|---|---|
2 | 3.98277278 | 4.04665506 | 4.04718956 | 0.0005345 | 0.00275556 |
4 | 4.03068449 | 4.04714980 | 4.04718956 | 0.00003976 | 0.00017222 |
8 | 4.04303347 | 4.04718692 | 4.04718956 | 0.00000264 | 0.00001076 |
16 | 4.04614856 | 4.04718939 | 4.04718956 | 0.00000017 | 0.00000067 |
32 | 4.04692918 | 4.04718955 | 4.04718956 | 0.00000001 | 0.00000004 |
64 | 4.04712446 | 4.04718956 | 4.04718956 | 0 | 0 |
20384 | 4.04718956 |
Здесь - коэффициент измельчения шага . Исходная величина шага .
На иллюстрации черная сплошная линия - исходная формула, красная пунктирная - экстраполированная.
Как мы видим, разница между экстраполированными и неэкстраполированными результатами значительна. Уже при величине шага в мы можем найти значение интеграла с точностью , тогда как в исходной формуле нам для достижения такой точности пришлось бы задать величину шага .
Рекомендации программисту
Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии. Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца.
Заключение
Мы получили более точные результаты при меньшем количестве вычислений функции, чем в исходном методе. Избавившись от степени, составляющей ошибку, мы пришли к результату, который в противном случае был бы недостижим без значительного уменьшения размера шага. Таким образом, мы преобразовали весьма посредственный алгоритм вычисления интегралов в довольно эффективный.
Список литературы
- А.А.Самарский, А.В.Гулин. Численные методы М.: Наука, 1989.
- Fundamental Methods of Numerical Extrapolation With Applications, mit.edu