Метод простых итераций

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

(Различия между версиями)
Перейти к: навигация, поиск
 
(69 промежуточных версий не показаны.)
Строка 1: Строка 1:
== Постановка задачи ==
== Постановка задачи ==
Пусть есть функция <tex>y = f(x)</tex>.<br>
Пусть есть функция <tex>y = f(x)</tex>.<br>
-
Требуется найти корень этой функции, то есть <tex>x</tex> при котором <tex>f(x)=0</tex><br>
+
Требуется найти корень этой функции: такой <tex>x</tex> при котором <tex>f(x)=0</tex><br>
-
Решение необходимо найти численно, то есть для реализации на ЭВМ. Для решения этой задачи предлагается использовать метод простых итераций и его обобщения.
+
Решение необходимо найти численно, то есть для реализации на ЭВМ. Для решения этой задачи предлагается использовать метод простых итераций.
-
=== Метод простых итераций в общем виде ===
+
== Метод простых итераций в общем виде ==
-
Заменеим исходное уравнение <tex>f(x)=0</tex> на эквивалентное <tex>g(x)=x</tex>.<br>
+
Заменим исходное уравнение <tex>f(x)=0</tex> на эквивалентное <tex>g(x)=x</tex>будем строить итерации по правилу <tex>x_{n+1} = g(x_n)</tex>. Таким образом метод простой итерации - это одношаговый итерационный процесс. Для того, что бы начать данный процесс, необходимо знать начальное приближение <tex>x_0</tex>. Выясним условия сходимости метода и выбор начального приближения.
-
Итерации будем строить по правилу <tex>g(x_n)=x_{n+1}</tex><br>
+
===Сходимость метода простых итераций===
-
Для сходимости метода очень важен выбор функции <tex>g(x)</tex>, поэтому ее обычно берут вида <tex>g(x)=x+s(x)f(x).</tex><br>
+
Метод сходится, если при <tex>k \to \infty </tex> последовательность {<tex>x_n</tex>} имеет предел.<br>
 +
Обозначим <tex>U_r(a)</tex> окресность точки <tex>a</tex> радиуса <tex>r</tex>, то есть <tex>U_r(a) = \{x:|x-a|<r\}</tex>.<br>
 +
'''Теорема 1.''' Если <tex>g(x)</tex> липшиц-непрерывна с константой <tex>q \in (0,1)</tex> на <tex>U_r(a)</tex>, то есть выполняется
 +
<center><tex>|g(x'')-g(x')|<q|x''-x'|</tex>,</center>
 +
при этом если также выполнено
 +
<center><tex>|g(a)-a|<(1-q)r</tex>,</center>то уравнение <tex>x = g(x)</tex> имеет единственное решение на <tex>U_r(a)</tex> и метод простой итерации сходится к решению при любом выборе начального приближения <tex>x_0 \in U_r(a)</tex>.Так же справедлива оценка:
 +
<center><tex>|x_k-x_*|<q^k|x_0-x_*|</tex>,</center>
 +
где <tex>x_*</tex> - точное решение.<br><br>
 +
Из оценки видно, что метод линеен.
 +
Пусть <tex>g(x)</tex> непрерывно дифференцируема на <tex>U_r(a)</tex>, тогда из теоремы вытекают следующие утверждения:<br>
 +
'''Следствие 1.''' Если <tex>|g'(x)| \le q < 1</tex> для <tex>x \in U_r(a)</tex>, выполнено <tex>|g(a)-a|<(1-q)r</tex>, и <tex>x_0 \in U_r(a)</tex>, тогда уравнение <tex>x = g(x)</tex> имеет единственное решение на <tex>U_r(a)</tex> и метод простой итерации сходится к решению.<br><br>
 +
'''Следствие 2.''' Если уравнение <tex>x = g(x)</tex> имеет решение <tex>x_*</tex>, <tex>g(x)</tex> непрерывно дифференцируема на <tex>U_r(x_*)</tex> и <tex>|g'(x_*)|<1</tex>. Тогда существует <tex>\eps > 0</tex> такое, что на <tex>U_{\eps}(x_*)</tex> уравнение не имеет других решений и метод простой итерации сходится к решению при <tex>x_0 \in U_{\eps}(x_*)</tex><br>
 +
===Геометрическая интерпретация===
 +
Рассмотрим график функции <tex> y = g(x)</tex>. Это озночает, что решение уравнения <tex>f(x) = 0</tex> и <tex>x=g(x)</tex> - это точка пересечения <tex>g(x)</tex> с прямой <tex>y = x</tex>:<br>
 +
<center>[[Изображение:PowerIterationMethod.jpg|300px]]</center>
 +
И следующая итерация <tex>x_{x+1} = g(x_n)</tex> - это координата <tex>x</tex> пересечения горизонтальной прямой точки <tex>(x_n g(x_n))</tex> с прямой <tex>y = x</tex>.<br>
 +
<center>[[Изображение:PowerIterationMethod4.jpg|300px]]</center>
 +
<br>
 +
Из рисунка наглядно видно требование сходимости <tex>|g'(x)|<1</tex>. Чем ближе производная <tex>g'(x)</tex> к <tex>0</tex>, тем быстрее сходится алгоритм. В зависимости от знака производной вблизи решения приближения могут строится по разному. Если <tex>g'(x)<0</tex>, то каждое следующее приближение строится с другой стороны от корня:<br>
 +
<center>[[Изображение:PowerIterationMethod2.jpg|300px]]</center>
 +
==Метод релаксации==
 +
Так как для сходимости метода очень важен выбор функции <tex>g(x)</tex>, ее обычно берут вида <center><tex>g(x)=x+s(x)f(x)</tex> <tex>(1)</tex>.</center>
Где <tex>s(x)</tex> не меняет знака на отрезке, на котором ищется корень функции.<br>
Где <tex>s(x)</tex> не меняет знака на отрезке, на котором ищется корень функции.<br>
-
Таким образом метод простой итерации - это одношаговый итерационны процесс.<br>
+
Положим <tex>s(x) = c = const </tex> и рассмотрим метод в этом случае.<br>
-
===Метод релаксации===
+
Тогда получим метод 'релаксации':
-
Положим <tex>s(x) = c = const </tex> b и рассмотрим метод в этом случае.<br>
+
<center><tex>f(x_n) = \frac{x_{n+1}-x_{n}}{c}</tex>,</center>
-
Тогда получим <tex>f(x_n) = \frac{x_{n+1}-x{n}}{c}</tex>
+
для которого <tex>g'(x) = 1+cf'(x)</tex>, и метод сходится при условии
-
===Сходимость метода===
+
<center><tex>-2<cf'(x_*)<0</tex>,</center>
 +
Пусть в некоторой окресности корня выполняются условия
 +
<center><tex>f'(x)<0, 0<m_1<|f'(x)|<M_1</tex>,</center>
 +
Тогда метод релаксации сходится при <tex>c \in (0,\frac{2}{M_1}).</tex>
 +
===Выбор параметра===
 +
Оценим погрешность метода релаксации <tex>z_k=x_k-x_*</tex>
 +
<center><tex>f(x_*+z_n) = \frac{z_{n+1}-z_{n}}{c}</tex>,</center>
 +
Применяя теорему о среднем получаем
 +
<center><tex>f'(x_*+{\theta}z_n)z_n = \frac{z_{n+1}-z_{n}}{c}</tex>,</center>
 +
Отсюда
 +
<center><tex>|z_{n+1}|\le|1+cf'(x_*+{\theta}z_n)||z_n|\le max|1+cf'(x_*+oz_n)||z_n|</tex>,</center>
 +
Следовательно
 +
<center><tex>|z_{n+1}|\le max\{|1-cM_1|,|1-cm_1|\}|z_n|</tex>,</center>
 +
Таким образом задача сводится к нахождению минимума функции <tex>F(c)</tex>
 +
<center><tex>F(c) = max\{|1-cM_1|,|1-cm_1|\}</tex>,</center>
 +
Из рассмотрения графика функции видно, что точка минимума определяется
 +
<center><tex>|1-cM_1| = |1-cm_1|</tex> при <tex>c\ne 0 </tex>,</center>
 +
и равна
 +
<center><tex>c = \frac{2}{M_1+m_1}</tex> </center>
 +
===Ускорение сходимости===
 +
Как следует из '''Теоремы 1''', метод простых итераций линеен, то есть
 +
<center><tex>x_n-x_* \approx aq^n</tex> </center>
 +
Воспользуемся этим для оценки погрешности на каждой итерации. Запомним 3 последние итерации и выпишем их оценки:
 +
<center><tex> x_{k}-x_{*} \approx aq^{k}</tex> </center>
 +
<center><tex> x_{k+1}-x_{*} \approx aq^{k+1}</tex> </center>
 +
<center><tex> x_{k+2}-x_{*} \approx aq^{k+2}</tex> </center>
 +
Где <tex> x_{k},x_{k+1},x_{k+2}</tex> нам известны (вычисленны по какому то линейному алгоритму),а <tex>a,q,x_*</tex> найдем из системы. Получим:
 +
<center><tex> x_{*} \approx x_{k+2} - \frac{(x_{k+2}-x_{k+1})^2}{x_{k+2}-2x_{k+1}+x_{k}}</tex> <tex>(2)</tex></center>
 +
Метод ускорения сходимости заключается в том, что после вычисления 3 приближений по линейно сходящемуся алгоритму, вычисляется новое приближение по уточняющему правилу (2).<br>
 +
Применительно к методу релаксации имеем:
 +
<center><tex> x_{*} \approx x_{k+2} - \frac{(x_{k+2}-x_{k+1})^2}{x_{k+2}-2x_{k+1}+x_{k}}</tex> </center><br>
 +
<center><tex>\frac{x_{n+1}-x_{n}}{c} = f(x_n)</tex> </center>
 +
Следовательно
 +
<center><tex> x_{*} \approx x_{k} - c\frac{f^2(x_k)}{f(x_k)-f(x_k-cf(x_k))}</tex> </center>
 +
Можно показать, что данный метод имеет уже квадратичную скорость сходимости.
 +
==Метод Вегстейна==
 +
Метод Вегстейна, вообще говоря, является модификацией метода секущих, однако его можно назвать и улучшенным методом простой итерации, преобразовав вычислительню формулу
 +
<center><tex>x_k = x_{k-1} - \frac{f(x_{k-1})(x_{k-1}-x_{k-2})}{f(x_{k-1})-f(x_{k-2})}</tex></center>
 +
к виду
 +
<center><tex>x_{k}=x_{k-1} - \frac{x_{k-1}-g(x_{k-1})}{(1-\frac{g(x_{k-1})-g(x_{k-2)}}{x_{k-1}-x_{k-2}})}</tex></center>
 +
Это двухшаговый метод, и для начала вычислений необходимо задать 2 приближения <tex>x_0,x_1</tex>.
 +
== Программная реализация ==
 +
Все методы были реализованы на языке C++. Доступ к методам осуществяется через класс
 +
PowerIterationMethod
 +
пример кода:
 +
PowerIterationMethod::PowerIterationParams *params =
 +
new PowerIterationMethod::PowerIterationParams (
 +
f1 // Исходная функция
 +
,s1 // Функция s(x) в формусле (1) или константа в методе релаксации
 +
,1 // Начальное приближение
 +
,0 // Второе приближение для метода Вегстейна
 +
,0 // Допустимая погрешность решения
 +
,1000 // Максимальное количество итераций
 +
);
 +
PowerIterationMethod *method = new PowerIterationMethod (params);
 +
method->simpleIteration (); // Вычисление по методу простой итерации
 +
printf ("%f\n",method->getResult ());
 +
printf ("%f",method->getEps ());
-
== Числовые примеры ==
+
== Примеры тестирования ==
-
== Рекомендации программисту ==
+
Ошибкой будем считать <tex>\eps = f(x_k)-f(x_*) = f(x_k)</tex> и проверим скорость сходимости методов относительно друг друга.<br>
 +
<tex>f(x) = -(\log(x)-\cos(x))</tex>
 +
<tex>\eps = 10^{-5}</tex><br>
 +
Начальное приближение <tex> x_0 = 1</tex><br>
 +
1. Метод простой итерации с <tex>s(x) = 1</tex>.<br>
 +
Сходимость за 28 шагов.<br>
 +
2. Метод простой итерации с <tex>s(x) = \sin x</tex>.<br>
 +
Сходимость за 21 шаг.<br>
 +
3. Ускоренный метод простой итерации.<br>
 +
Сходимость за 3 шага.<br>
 +
4. Метод Вегстейна.<br>
 +
Сходимость за 3 шага.<br><br><br>
 +
<tex>f(x) = -(\frac{1}{2}+x^2-\cos x)</tex>
 +
Корень <tex>x_* \approx 0.58 </tex><br>
 +
<tex>\eps = 10^{-5}</tex><br>
 +
Начальное приближение <tex> x_0 = 1</tex><br>
 +
1. Метод простой итерации с <tex>s(x) = 1</tex>.<br>
 +
Сходимость за 23 шагов.<br>
 +
2. Метод простой итерации с <tex>s(x) = \sin x</tex>.<br>
 +
Сходимость за 5 шаг.<br>
 +
3. Ускоренный метод простой итерации.<br>
 +
Сходимость за 4 шага.<br>
 +
4. Метод Вегстейна.<br>
 +
Сходимость за 4 шага.<br><br><br>
 +
<tex>f(x) = -(\frac{1}{2}+x^2-\cos x)</tex>
 +
Корень <tex>x_* \approx 0.58 </tex><br>
 +
<tex>\eps = 10^{-8}</tex><br>
 +
Начальное приближение <tex> x_0 = 0.4</tex><br>
 +
1. Метод простой итерации с <tex>s(x) = 1</tex>.<br>
 +
Сходимость за 43 шагов.<br>
 +
2. Метод простой итерации с <tex>s(x) = \sin x</tex>.<br>
 +
Сходимость за 7 шагов.<br>
 +
3. Ускоренный метод простой итерации.<br>
 +
Сходимость за 5 шагов.<br>
 +
4. Метод Вегстейна.<br>
 +
Сходимость за 7 шагов.<br>
 +
 
 +
 
 +
Исходный код можно скачать [[Media:PowerIterationMethod.zip|Код программы]]
== Заключение ==
== Заключение ==
== Ссылки ==
== Ссылки ==
* [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
* [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
-
 
== Список литературы ==
== Список литературы ==
-
* ''А.А.Самарский, А.В.Гулин.''&nbsp; Численные методы. Москва «Наука», 1989.
+
* ''А.А.Самарский, А.В.Гулин.'' Численные методы. Москва «Наука», 1989.
-
* ''Н.Н.Калиткин.''&nbsp; Численные методы. Москва «Наука», 1978.
+
* ''Н.Н.Калиткин.'' Численные методы. Москва «Наука», 1978.
{{stub}}
{{stub}}

Текущая версия

Содержание

Постановка задачи

Пусть есть функция y = f(x).
Требуется найти корень этой функции: такой x при котором f(x)=0
Решение необходимо найти численно, то есть для реализации на ЭВМ. Для решения этой задачи предлагается использовать метод простых итераций.

Метод простых итераций в общем виде

Заменим исходное уравнение f(x)=0 на эквивалентное g(x)=x,и будем строить итерации по правилу x_{n+1} = g(x_n). Таким образом метод простой итерации - это одношаговый итерационный процесс. Для того, что бы начать данный процесс, необходимо знать начальное приближение x_0. Выясним условия сходимости метода и выбор начального приближения.

Сходимость метода простых итераций

Метод сходится, если при k \to \infty последовательность {x_n} имеет предел.
Обозначим U_r(a) окресность точки a радиуса r, то есть U_r(a) = \{x:|x-a|<r\}.
Теорема 1. Если g(x) липшиц-непрерывна с константой q \in (0,1) на U_r(a), то есть выполняется

|g(x'')-g(x')|<q|x''-x'|,

при этом если также выполнено

|g(a)-a|<(1-q)r,
то уравнение x = g(x) имеет единственное решение на U_r(a) и метод простой итерации сходится к решению при любом выборе начального приближения x_0 \in U_r(a).Так же справедлива оценка:
|x_k-x_*|<q^k|x_0-x_*|,

где x_* - точное решение.

Из оценки видно, что метод линеен. Пусть g(x) непрерывно дифференцируема на U_r(a), тогда из теоремы вытекают следующие утверждения:
Следствие 1. Если |g'(x)| \le q < 1 для x \in U_r(a), выполнено |g(a)-a|<(1-q)r, и x_0 \in U_r(a), тогда уравнение x = g(x) имеет единственное решение на U_r(a) и метод простой итерации сходится к решению.

Следствие 2. Если уравнение x = g(x) имеет решение x_*, g(x) непрерывно дифференцируема на U_r(x_*) и |g'(x_*)|<1. Тогда существует \eps > 0 такое, что на U_{\eps}(x_*) уравнение не имеет других решений и метод простой итерации сходится к решению при x_0 \in U_{\eps}(x_*)

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

Рассмотрим график функции  y = g(x). Это озночает, что решение уравнения f(x) = 0 и x=g(x) - это точка пересечения g(x) с прямой y = x:

И следующая итерация x_{x+1} = g(x_n) - это координата x пересечения горизонтальной прямой точки (x_n g(x_n)) с прямой y = x.


Из рисунка наглядно видно требование сходимости |g'(x)|<1. Чем ближе производная g'(x) к 0, тем быстрее сходится алгоритм. В зависимости от знака производной вблизи решения приближения могут строится по разному. Если g'(x)<0, то каждое следующее приближение строится с другой стороны от корня:

Метод релаксации

Так как для сходимости метода очень важен выбор функции g(x), ее обычно берут вида
g(x)=x+s(x)f(x) (1).

Где s(x) не меняет знака на отрезке, на котором ищется корень функции.
Положим s(x) = c = const и рассмотрим метод в этом случае.
Тогда получим метод 'релаксации':

f(x_n) = \frac{x_{n+1}-x_{n}}{c},

для которого g'(x) = 1+cf'(x), и метод сходится при условии

-2<cf'(x_*)<0,

Пусть в некоторой окресности корня выполняются условия

f'(x)<0, 0<m_1<|f'(x)|<M_1,

Тогда метод релаксации сходится при c \in (0,\frac{2}{M_1}).

Выбор параметра

Оценим погрешность метода релаксации z_k=x_k-x_*

f(x_*+z_n) = \frac{z_{n+1}-z_{n}}{c},

Применяя теорему о среднем получаем

f'(x_*+{\theta}z_n)z_n = \frac{z_{n+1}-z_{n}}{c},

Отсюда

|z_{n+1}|\le|1+cf'(x_*+{\theta}z_n)||z_n|\le max|1+cf'(x_*+oz_n)||z_n|,

Следовательно

|z_{n+1}|\le max\{|1-cM_1|,|1-cm_1|\}|z_n|,

Таким образом задача сводится к нахождению минимума функции F(c)

F(c) = max\{|1-cM_1|,|1-cm_1|\},

Из рассмотрения графика функции видно, что точка минимума определяется

|1-cM_1| = |1-cm_1| при c\ne 0 ,

и равна

c = \frac{2}{M_1+m_1}

Ускорение сходимости

Как следует из Теоремы 1, метод простых итераций линеен, то есть

x_n-x_* \approx aq^n

Воспользуемся этим для оценки погрешности на каждой итерации. Запомним 3 последние итерации и выпишем их оценки:

 x_{k}-x_{*} \approx aq^{k}
 x_{k+1}-x_{*} \approx aq^{k+1}
 x_{k+2}-x_{*} \approx aq^{k+2}

Где  x_{k},x_{k+1},x_{k+2} нам известны (вычисленны по какому то линейному алгоритму),а a,q,x_* найдем из системы. Получим:

 x_{*} \approx x_{k+2} - \frac{(x_{k+2}-x_{k+1})^2}{x_{k+2}-2x_{k+1}+x_{k}} (2)

Метод ускорения сходимости заключается в том, что после вычисления 3 приближений по линейно сходящемуся алгоритму, вычисляется новое приближение по уточняющему правилу (2).
Применительно к методу релаксации имеем:

 x_{*} \approx x_{k+2} - \frac{(x_{k+2}-x_{k+1})^2}{x_{k+2}-2x_{k+1}+x_{k}}

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

Следовательно

 x_{*} \approx x_{k} - c\frac{f^2(x_k)}{f(x_k)-f(x_k-cf(x_k))}

Можно показать, что данный метод имеет уже квадратичную скорость сходимости.

Метод Вегстейна

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

x_k = x_{k-1} - \frac{f(x_{k-1})(x_{k-1}-x_{k-2})}{f(x_{k-1})-f(x_{k-2})}

к виду

x_{k}=x_{k-1} - \frac{x_{k-1}-g(x_{k-1})}{(1-\frac{g(x_{k-1})-g(x_{k-2)}}{x_{k-1}-x_{k-2}})}

Это двухшаговый метод, и для начала вычислений необходимо задать 2 приближения x_0,x_1.

Программная реализация

Все методы были реализованы на языке C++. Доступ к методам осуществяется через класс

PowerIterationMethod

пример кода:

PowerIterationMethod::PowerIterationParams *params = 
        new PowerIterationMethod::PowerIterationParams (
        f1   // Исходная функция
       ,s1   // Функция s(x) в формусле (1) или константа в методе релаксации
       ,1    // Начальное приближение
       ,0    // Второе приближение для метода Вегстейна
       ,0    // Допустимая погрешность решения
       ,1000 // Максимальное количество итераций
       );
PowerIterationMethod *method = new PowerIterationMethod (params);
method->simpleIteration (); // Вычисление по методу простой итерации
printf ("%f\n",method->getResult ());
printf ("%f",method->getEps ());

Примеры тестирования

Ошибкой будем считать \eps = f(x_k)-f(x_*) = f(x_k) и проверим скорость сходимости методов относительно друг друга.
f(x) = -(\log(x)-\cos(x)) \eps = 10^{-5}
Начальное приближение  x_0 = 1
1. Метод простой итерации с s(x) = 1.
Сходимость за 28 шагов.
2. Метод простой итерации с s(x) = \sin x.
Сходимость за 21 шаг.
3. Ускоренный метод простой итерации.
Сходимость за 3 шага.
4. Метод Вегстейна.
Сходимость за 3 шага.


f(x) = -(\frac{1}{2}+x^2-\cos x) Корень x_* \approx 0.58
\eps = 10^{-5}
Начальное приближение  x_0 = 1
1. Метод простой итерации с s(x) = 1.
Сходимость за 23 шагов.
2. Метод простой итерации с s(x) = \sin x.
Сходимость за 5 шаг.
3. Ускоренный метод простой итерации.
Сходимость за 4 шага.
4. Метод Вегстейна.
Сходимость за 4 шага.


f(x) = -(\frac{1}{2}+x^2-\cos x) Корень x_* \approx 0.58
\eps = 10^{-8}
Начальное приближение  x_0 = 0.4
1. Метод простой итерации с s(x) = 1.
Сходимость за 43 шагов.
2. Метод простой итерации с s(x) = \sin x.
Сходимость за 7 шагов.
3. Ускоренный метод простой итерации.
Сходимость за 5 шагов.
4. Метод Вегстейна.
Сходимость за 7 шагов.


Исходный код можно скачать Код программы

Заключение

Ссылки

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

  • А.А.Самарский, А.В.Гулин. Численные методы. Москва «Наука», 1989.
  • Н.Н.Калиткин. Численные методы. Москва «Наука», 1978.
Личные инструменты