Участник:Gukov/Песочница

Материал из MachineLearning.

Перейти к: навигация, поиск

Содержание

Введение

Метод Ньютона-Рафсона

Пусть имеется некоторая функция F(x) и необходимо найти такие значения аргумента x, для которых

(1)
F(x)=0

Перепишем (1) в виде

(1.1)
x = f(x)

и запишем n+1-ое приближения корня (1.1), при этом делая поправку \alpha к очередному значению x_n

(2)
x_{n+1} = x_n + \alpha \Delta x_n,

где \Delta x_n = f(x_n) - x_n, и положим \alpha = \frac{1}{1-f'(x_n)}.

Тогда (2) перепишется в виде

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

Нетрудно видеть, что (3) эквивалентно простому методу последовательных приближений

x_{n+1} = g(x_n),

где g(x) = \frac{f(x) - x f'(x)}{1 - f'(x)}.

Вспомним также, что если |g'(x)| < 1, то метод сходится. Имеем

g'(x) = \frac{f''(x)[f(x) - x]}{[1 - f'(x)]^2}.

Но так как справедливо соотношение (1.1), то для x, достаточно близких к решению (1.1), выражение в скобках в числителе дроби становится малым. Поэтому итерационный метод, описываемый формулой (3) сходится, если

  • 1. Начальное приближение x_0 выбрано достаточно близко к решению x = f(x).
  • 2. Производная f''(x) не становится слишком большой.
  • 3. Производная f'(x) не слишком близка к 1.

Это и есть метод Ньютона-Рафсона. Обычно его записывают в виде

(4)
x_{n+1} = x_{n} - \frac{F(x_n)}{F'(x_n)},

где F(x) = f(x) - x = 0

Таким образом, мы вернулись к уравнению в форме (1), и условия сходимости принимают следующий вид

  • 1. Начальное приближение x_0 выбрано достаточно близко к корню уравнения F(x) = 0.
  • 2. Производная F''(x) не становится очень большой.
  • 3. Производная F'(x) не слишком близка к 0.

Случай почти равных корней

Рисунок 1
Рисунок 1

Условие 3 сходимости метода Ньютона-Рафсона означает, что никакие два корня не находятся слишком близко один к другому. Соответствующая ситуация представлена на рисунке 1 (масштаб сильно увеличен). Заметим, что производная f'(x) близка к 1 при x, равном обоим значениям корней, a_1 и a_2. Более того, на основании теоремы Лагранжа, можно утверждать, что f'(x) = 1 где-то между a_1 и a_2.

Рассмотрим, что случится, если принять x_0 в качестве исходного значения для корня a_1. Касательная, проведенная через точку C, пересечет прямую y=x в точке A, и следующее приближение будет равно x_1. Касательная, проведенная через точку B, пересекает прямую в точке D, и в качестве следующего приближеня получается снова x_0. Итерационный процесс, таким образом, осциллирует между x_0 и x_1 до бесконечности, не сходясь ни к одному значению корня. Иначе говоря, не удается отделить эти два корня, потому что они расположены слишком близко.

Изложение метода

Экстраполяция Ричардсона

Предположим, что для вычисления интеграла (1) отрезок [a, b] разбит на N равных отрезков длины h = \frac{b-a}N и на каждом частичном отрезке применяется одна и та жа квадратурная формула. Тогда исходный интеграл I заменяется некоторой квадратурной суммой I_h, причем возникающая погрешность зависит от шага сетки h. Для некоторых квадратурных формул удается получить разложение погрешности I_h - I по степеням h. Предположим, что для данной квадратурной суммы I_h существует разложение:

( 3)
I_h = I_0 + a_1 h^{\alpha _1} + a_2 h^{\alpha _2} + \ldots + a_m h^{\alpha _m} + O(h^{\alpha _{m+1}}),

где 0 < \alpha _1 < \alpha _2 < \ldots < \alpha _m < \alpha _{m+1} и коэффициенты \{ a_i \} \subset \mathbb{R} не зависят от h. При этом величины \{ \alpha _i \} \subset \mathbb{R} предполагаются известными. Теперь предположим:

I(\frac{h}{r}) = I_0 + a_1\,\frac{h^{\alpha _1}}{r^{\alpha 1}} + a_2\,\frac{h^{\alpha _2}}{r^{\alpha 2}} + \ldots

Чтобы избавиться от степени h^{\alpha _1}, составляющей ошибку (ибо среди всех слагаемых, составляющих ошибку, слагамое при h^{\alpha _1} является наибольшим) вычислим величину r^{\alpha _1}\,I(\frac{h}r) - I(h). Имеем:

r^{\alpha _1}\,I(\frac{h}{r}) - I(h) = r^{\alpha _1}\,I_0 - I_0 + \fbox{a_1 h^{\alpha _1} - a_1 h^{\alpha _1}} + a_2\,\frac{h^{\alpha _2}}{r^{\alpha _2 - \alpha _1}} - a_2 h^{\alpha _2} + \ldots

Отсюда

I_1(h) = \frac{r^{\alpha _1}I(\frac{h}r) - I(h)}{r^{\alpha _1} - 1} = I_0 + \underbrace{a_2\,\frac{r^{\alpha _1 - \alpha _2} - 1}{r^{\alpha _1} - 1}}_{a_2 '}\,h^{\alpha _2} + \ldots

то есть имеем более точное приближение к интегралу I.

Таким образом, рекуррентную формулу можно записать в виде:

( 4)
I_{i+1}(h) = I_i(\frac{h}r) + \frac{I_i(\frac{h}r) - I_i(h)}{r^{\alpha _1} - 1}

Заметим, что r - величина, на которую мы делим размер шага при каждом новом вычислении I. Разумно положить s = 2, т.к. большие значения s могут вызвать резкое увеличение количества вычислений.

Для наглядности представим процесс экстраполирования следующей таблицей:


\line(1, 0){200}

I_0^{0}(h)
I_0^{1}(\frac{h}r) I_1^{0}(h)
I_0^{2}(\frac{h}{r^2}) I_1^{1}(\frac{h}r) I_2^{0}(h)
\vdots \vdots \vdots \ddots

\line(1,0){200}

Метод Рунге

Пусть существует асимптотическое разложение вида:

D_N(f) = I_N(f) - I_0 = \alpha _2 h^2 + \alpha _4 h^4 + \alpha _6 h^6 + \ldots,

где f(x) - достаточно гладкая функция, а hN = b - a.

Проведем расчеты на двух равномерных сетках с шагами h_1 и h_2 соответственно и найдем выражения I^{h_1}[f] = I_{N_1}[f] и I^{h_2}[f] = I_{N_2}[f], h_1 N_1 = h_2 N_2 = b - a. Потребуем, чтобы погрешность для их линейной комбинации:

{\tilde D}^{h}(f) = \sigma D^{h_1}(f) + (1 - \sigma)D^{h_2}(f)

была величиной более высокого порядка по сравнению с D^{h_1} и D^{h_2}. Если для D^h = D_N имеет место формула вида

 D^{h} = I^{h}[f] - I[f] = \alpha _p h^p + \alpha _q h^q + \ldots,     q > p,

то для

{\tilde D}^{h} = (\sigma I^{h_1}[f] + (1 - \sigma)I^{h_2}[f]) - I[f]

получим

{\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

Выберем параметр \sigma из условия \sigma h_1^{p} + (1 - \sigma)h_2^{p}:

\sigma = \frac{h_2^p}{h_2^p - h_1^p}.

Имеем

{\tilde D}^h(f) = \alpha _q(\sigma h_1^q + (1 - \sigma)h_2^q) + \ldots = O(h^q),       h = max(h_1, h_2),

причем \sigma h_1^q + (1 - \sigma)h^q_2 < 0. Так, если p = 2, q = 4, то {\tilde D}^{h}(f) = - \alpha _4 h_1^2 h_2^2 + \ldots = O(h^4). Таким образом, проведя вычисления на двух сетках с шагами h_1 и h_2 \ne h_1, мы повысили порядок точности на 2 (на (q - p)) для {\tilde I} = \sigma I^{h_1} + (1 - \sigma)I^{h_2}.

Процесс Эйткена

Метод расчета на нескольких сетках применяется для повышения порядка точности даже в том случае, когда неизвестен порядок главного члена погрешности. Предположим, чт для погрешности имеет место представление

D^{h}(f) = \alpha _p h^{p} + \alpha _q h^{q} + \ldots,     q > p,

так что

I^{h}[f] = I[f] + \alpha _p h^p + \alpha _q h^{q} + \ldots

Проведем вычисления на трех сетках: h_1 = h, h_2 = \rho h, h_3 = \rho^{2} h (0 < \rho < 1). Определим p. При этом пренебрегаем членом O(h^q). Образуем отношение

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

и найдем

p \approx \frac{\ln A}{\ln(\frac{1}{\rho})}.

Зная приближенное значение p, можно методом методом Рунге повысить порядок точности. Для этого образуем комбинацию {\tilde I}^h = \sigma I^{h_1} + (1 - \sigma)I^{h_2} и выберем \sigma так, чтобы  \sigma h_1^p+(1 - \sigma)h_2^p = (\sigma + (1 - \sigma)\rho^p)h^p = 0 , то есть \sigma = \frac{\rho^p}{\rho^p - 1} = \frac{1}{1 - A}. Тогда для погрешности {\tilde D}^h = {\tilde I}^h - I получаем

 {\tilde D}^h(f) = O(h^q) .

Числовой пример

Найдем с помощью квадратурной формулы трапеций приближенное значение интеграла, применив экстраполяцию Ричардсона (данный метод называется методом Ромберга):

\int^{5}_{1} \ln x\, dx = x \ln x - x |_1^{5} = 5\ln 5 - 4

В нижеследующей таблице представлены результаты работы программы:

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

Здесь r - коэффициент измельчения шага h. Исходная величина шага h = 2.

На иллюстрации черная сплошная линия - исходная формула, красная пунктирная - экстраполированная.

Как мы видим, разница между экстраполированными и неэкстраполированными результатами значительна. Уже при величине шага в h = \frac{1}{64} мы можем найти значение интеграла с точностью 10^{-8}, тогда как в исходной формуле нам для достижения такой точности пришлось бы задать величину шага h = \frac{1}{10192}.

Рекомендации программисту

Программируя экстраполяцию Ричардсона следует предпочесть итерацию рекурсии. Также реализуя многократное экстраполирование итеративно, нам не нужно хранить всю таблицу значений - достаточно иметь в распоряжении два последних столбца.

Заключение

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

Список литературы

Личные инструменты