Метод простых итераций
Материал из MachineLearning.
(44 промежуточные версии не показаны) | |||
Строка 1: | Строка 1: | ||
== Постановка задачи == | == Постановка задачи == | ||
Пусть есть функция <tex>y = f(x)</tex>.<br> | Пусть есть функция <tex>y = f(x)</tex>.<br> | ||
- | Требуется найти корень этой функции | + | Требуется найти корень этой функции: такой <tex>x</tex> при котором <tex>f(x)=0</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>k \to \infty </tex> последовательность {<tex>x_n</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> | Обозначим <tex>U_r(a)</tex> окресность точки <tex>a</tex> радиуса <tex>r</tex>, то есть <tex>U_r(a) = \{x:|x-a|<r\}</tex>.<br> | ||
- | '''Теорема.''' Если <tex>g(x)</tex> липшиц-непрерывна с константой <tex>q \in (0,1)</tex> на <tex>U_r(a)</tex>, то есть выполняется | + | '''Теорема 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(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> | + | <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> | <center><tex>|x_k-x_*|<q^k|x_0-x_*|</tex>,</center> | ||
где <tex>x_*</tex> - точное решение.<br><br> | где <tex>x_*</tex> - точное решение.<br><br> | ||
+ | Из оценки видно, что метод линеен. | ||
Пусть <tex>g(x)</tex> непрерывно дифференцируема на <tex>U_r(a)</tex>, тогда из теоремы вытекают следующие утверждения:<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> | '''Следствие 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> | '''Следствие 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>, ее обычно берут вида <tex>g(x)=x+s(x)f(x)</tex>. Где <tex>s(x)</tex> не меняет знака на отрезке, на котором ищется корень функции.<br> | + | Так как для сходимости метода очень важен выбор функции <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) = c = const </tex> и рассмотрим метод в этом случае.<br> | Положим <tex>s(x) = c = const </tex> и рассмотрим метод в этом случае.<br> | ||
- | Тогда получим <tex>f(x_n) = \frac{x_{n+1}-x_{n}}{c}</tex>. <br> | + | Тогда получим метод 'релаксации': |
+ | <center><tex>f(x_n) = \frac{x_{n+1}-x_{n}}{c}</tex>,</center> | ||
+ | для которого <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]] | ||
- | |||
== Список литературы == | == Список литературы == | ||
* ''А.А.Самарский, А.В.Гулин.'' Численные методы. Москва «Наука», 1989. | * ''А.А.Самарский, А.В.Гулин.'' Численные методы. Москва «Наука», 1989. |
Текущая версия
Содержание |
Постановка задачи
Пусть есть функция .
Требуется найти корень этой функции: такой при котором
Решение необходимо найти численно, то есть для реализации на ЭВМ. Для решения этой задачи предлагается использовать метод простых итераций.
Метод простых итераций в общем виде
Заменим исходное уравнение на эквивалентное ,и будем строить итерации по правилу . Таким образом метод простой итерации - это одношаговый итерационный процесс. Для того, что бы начать данный процесс, необходимо знать начальное приближение . Выясним условия сходимости метода и выбор начального приближения.
Сходимость метода простых итераций
Метод сходится, если при последовательность {} имеет предел.
Обозначим окресность точки радиуса , то есть .
Теорема 1. Если липшиц-непрерывна с константой на , то есть выполняется
при этом если также выполнено
где - точное решение.
Из оценки видно, что метод линеен.
Пусть непрерывно дифференцируема на , тогда из теоремы вытекают следующие утверждения:
Следствие 1. Если для , выполнено , и , тогда уравнение имеет единственное решение на и метод простой итерации сходится к решению.
Следствие 2. Если уравнение имеет решение , непрерывно дифференцируема на и . Тогда существует такое, что на уравнение не имеет других решений и метод простой итерации сходится к решению при
Геометрическая интерпретация
Рассмотрим график функции . Это озночает, что решение уравнения и - это точка пересечения с прямой :
И следующая итерация - это координата пересечения горизонтальной прямой точки с прямой .
Из рисунка наглядно видно требование сходимости . Чем ближе производная к , тем быстрее сходится алгоритм. В зависимости от знака производной вблизи решения приближения могут строится по разному. Если , то каждое следующее приближение строится с другой стороны от корня:
Метод релаксации
Так как для сходимости метода очень важен выбор функции , ее обычно берут видаГде не меняет знака на отрезке, на котором ищется корень функции.
Положим и рассмотрим метод в этом случае.
Тогда получим метод 'релаксации':
для которого , и метод сходится при условии
Пусть в некоторой окресности корня выполняются условия
Тогда метод релаксации сходится при
Выбор параметра
Оценим погрешность метода релаксации
Применяя теорему о среднем получаем
Отсюда
Следовательно
Таким образом задача сводится к нахождению минимума функции
Из рассмотрения графика функции видно, что точка минимума определяется
и равна
Ускорение сходимости
Как следует из Теоремы 1, метод простых итераций линеен, то есть
Воспользуемся этим для оценки погрешности на каждой итерации. Запомним 3 последние итерации и выпишем их оценки:
Где нам известны (вычисленны по какому то линейному алгоритму),а найдем из системы. Получим:
Метод ускорения сходимости заключается в том, что после вычисления 3 приближений по линейно сходящемуся алгоритму, вычисляется новое приближение по уточняющему правилу (2).
Применительно к методу релаксации имеем:
Следовательно
Можно показать, что данный метод имеет уже квадратичную скорость сходимости.
Метод Вегстейна
Метод Вегстейна, вообще говоря, является модификацией метода секущих, однако его можно назвать и улучшенным методом простой итерации, преобразовав вычислительню формулу
к виду
Это двухшаговый метод, и для начала вычислений необходимо задать 2 приближения .
Программная реализация
Все методы были реализованы на языке 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 ());
Примеры тестирования
Ошибкой будем считать и проверим скорость сходимости методов относительно друг друга.
Начальное приближение
1. Метод простой итерации с .
Сходимость за 28 шагов.
2. Метод простой итерации с .
Сходимость за 21 шаг.
3. Ускоренный метод простой итерации.
Сходимость за 3 шага.
4. Метод Вегстейна.
Сходимость за 3 шага.
Корень
Начальное приближение
1. Метод простой итерации с .
Сходимость за 23 шагов.
2. Метод простой итерации с .
Сходимость за 5 шаг.
3. Ускоренный метод простой итерации.
Сходимость за 4 шага.
4. Метод Вегстейна.
Сходимость за 4 шага.
Корень
Начальное приближение
1. Метод простой итерации с .
Сходимость за 43 шагов.
2. Метод простой итерации с .
Сходимость за 7 шагов.
3. Ускоренный метод простой итерации.
Сходимость за 5 шагов.
4. Метод Вегстейна.
Сходимость за 7 шагов.
Исходный код можно скачать Код программы
Заключение
Ссылки
Список литературы
- А.А.Самарский, А.В.Гулин. Численные методы. Москва «Наука», 1989.
- Н.Н.Калиткин. Численные методы. Москва «Наука», 1978.