Метод простой итерации, называемый также методом последовательного приближения, - это математический алгоритм нахождения значения неизвестной величины путем постепенного ее уточнения. Суть этого метода в том, что, как видно из названия, постепенно выражая из начального приближения последующие, получают все более уточненные результаты. Этот метод используется для поиска значения переменной в заданной функции, а также при решении систем уравнений, как линейных, так и нелинейных.
Рассмотрим, как данный метод реализуется при решении СЛАУ. Метод простой итерации имеет следующий алгоритм:
1. Проверка выполнения условия сходимости в исходной матрице. Теорема о сходимости: если исходная матрица системы имеет диагональное преобладание (т.е, в каждой строке элементы главной диагонали должны быть больше по модулю, чем сумма элементов побочных диагоналей по модулю), то метод простых итераций - сходящийся.
2. Матрица исходной системы не всегда имеет диагональное преобладание. В таких случаях систему можно преобразовать. Уравнения, удовлетворяющие условию сходимости, оставляют нетронутыми, а с неудовлетворяющими составляют линейные комбинации, т.е. умножают, вычитают, складывают уравнения между собой до получения нужного результата.
Если в полученной системе на главной диагонали находятся неудобные коэффициенты, то к обеим частям такого уравнения прибавляют слагаемые вида с i *x i, знаки которых должны совпадать со знаками диагональных элементов.
3. Преобразование полученной системы к нормальному виду:
x - =β - +α*x -
Это можно сделать множеством способов, например, так: из первого уравнения выразить х 1 через другие неизвестные, из второго- х 2 , из третьего- х 3 и т.д. При этом используем формулы:
α ij = -(a ij / a ii)
i = b i /a ii
Следует снова убедиться, что полученная система нормального вида соответствует условию сходимости:
∑ (j=1) |α ij |≤ 1, при этом i= 1,2,...n
4. Начинаем применять, собственно, сам метод последовательных приближений.
x (0) - начальное приближение, выразим через него х (1) , далее через х (1) выразим х (2) . Общая формула а матричном виде выглядит так:
x (n) = β - +α*x (n-1)
Вычисляем, пока не достигнем требуемой точности:
max |x i (k)-x i (k+1) ≤ ε
Итак, давайте разберем на практике метод простой итерации. Пример:
Решить СЛАУ:
4,5x1-1.7x2+3.5x3=2
3.1x1+2.3x2-1.1x3=1
1.8x1+2.5x2+4.7x3=4 с точностью ε=10 -3
Посмотрим, преобладают ли по модулю диагональные элементы.
Мы видим что условию сходимости удовлетворяет лишь третье уравнение. Первое и второе преобразуем, к первому уравнению прибавим второе:
7,6x1+0.6x2+2.4x3=3
Из третьего вычтем первое:
2,7x1+4.2x2+1.2x3=2
Мы преобразовали исходную систему в равноценную:
7,6x1+0.6x2+2.4x3=3
-2,7x1+4.2x2+1.2x3=2
1.8x1+2.5x2+4.7x3=4
Теперь приведем систему к нормальному виду:
x1=0.3947-0.0789x2-0.3158x3
x2=0.4762+0.6429x1-0.2857x3
x3= 0.8511-0.383x1-0.5319x2
Проверяем сходимость итерационного процесса:
0.0789+0.3158=0,3947 ≤ 1
0.6429+0.2857=0.9286 ≤ 1
0.383+ 0.5319= 0.9149 ≤ 1 , т.е. условие выполняется.
0,3947
Начальное приближение х (0) = 0,4762
0,8511
Подставляем данные значения в уравнение нормального вида, получаем следующие значения:
0,08835
x (1) = 0,486793
0,446639
Подставляем новые значения, получаем:
0,215243
x (2) = 0,405396
0,558336
Продолжаем вычисления до того момента, пока не приблизимся к значениям, удовлетворяющим заданному условию.
x (7) = 0,441091
Проверим правильность полученных результатов:
4,5*0,1880 -1.7*0,441+3.5*0,544=2,0003
3.1*0,1880+2.3*0,441-1.1x*0,544=0,9987
1.8*0,1880+2.5*0,441+4.7*0,544=3,9977
Результаты, полученные при подстановке найденных значений в исходные уравнения, полностью удовлетворяют условиям уравнения.
Как мы видим, метод простой итерации дает довольно точные результаты, однако для решения этого уравнения нам пришлось потратить много времени и проделать громоздкие вычисления.
Метод простых итераций основан на замене исходного уравнения эквивалентным уравнением:
Пусть известно начальное приближение к корню х = х 0 . Подставив его в правую часть уравнения (2.7), получим новое приближение , затем аналогичным образом получим и т. д.:
. (2.8)
Не при всех условиях итерационный процесс сходится к корню уравнения х . Рассмотрим этот процесс подробнее. На рис.2.6 приведена графическая интерпретация одностороннего сходящегося и расходящегося процесса. На рис.2.7 изображены двухсторонний сходящийся и расходящийся процессы. Расходящийся процесс характеризуется быстрым нарастанием значений аргумента и функции и аварийным завершением соответствующей программы.
При двухстороннем процессе возможно зацикливание, то есть бесконечное повторение одних и тех же значений функции и аргумента. Зацикливание отделяет расходящийся процесс от сходящегося.
Из графиков видно, что как при одностороннем, так и при двухстороннем процессе сходимость к корню определяется наклоном кривой вблизи корня. Чем меньше наклон, тем лучше сходимость. Как известно, тангенс угла наклона кривой равен производной кривой в данной точке.
Следовательно, чем меньше вблизи корня, тем быстрее сходится процесс.
Для того чтобы итерационный процесс был сходящимся, необходимо в окрестности корня выполнение следующего неравенства:
Переход от уравнения (2.1) к уравнению (2.7) можно осуществить различными способами в зависимости от вида функции f(x). При таком переходе необходимо построить функцию так, чтобы выполнялось условие сходимости (2.9).
Рассмотрим один из общих алгоритмов перехода от уравнения (2.1) к уравнению (2.7).
Умножим левую и правую части уравнения (2.1) на произвольную константу b и добавим к обеим частям неизвестное х. При этом корни исходного уравнения не изменятся:
Введем обозначение и перейдем от соотношения (2.10) к уравнению (2.8).
Произвольный выбор константы b позволит обеспечить выполнение условия сходимости (2.9). Критерием окончания итерационного процесса будет условие (2.2). На рис.2.8 приведена графическая интерпретация метода простых итераций при описанном способе представления (масштабы по осям X и Y различны).
Если функция выбрана в виде , то производная от этой функции будет . Наибольшая скорость сходимости будет при , тогда и итерационная формула (2.11) переходит в формулу Ньютона . Таким образом, метод Ньютона имеет самую высокую степень сходимости из всех итерационных процессов.
Программная реализация метода простых итераций выполнена в виде процедуры-подпрограммы Iteras (ПРОГРАММА 2.1).
Вся процедура практически состоит из одного цикла Repeat ... Until, реализующего формулу (2.11) с учетом условия прекращения итерационного процесса (формула (2.2)).
В процедуру встроена защита от зацикливания путем подсчета числа циклов с помощью переменной Niter. На практических занятиях необходимо убедиться путем прогона программы в том, как сказывается выбор коэффициента b и начального приближения на процессе поиска корня. При изменении коэффициента b характер итерационного процесса для исследуемой функции меняется. Он становится сначала двухсторонним, а потом зацикливается (рис.2.9). Масштабы по осям X и Y различны. Еще большее значение модуля b приводит к расходящемуся процессу.
Сравнение методов приближенного решения уравнений
Сравнение описанных выше методов численного решения уравнений проводилось с помощью программы, позволяющей на экране ПЭВМ наблюдать процесс нахождения корня в графическом виде. Процедуры, входящие в данную программу и реализующие сравниваемые методы, приведены ниже (ПРОГРАММА 2.1).
Рис. 2.3-2.5, 2.8, 2.9 являются копиями экрана ПЭВМ при окончании итерационного процесса.
В качестве исследуемой функции во всех случаях было взято квадратное уравнение x 2 -x-6 = 0, имеющее аналитическое решение х 1 = -2 и х 2 = 3. Погрешность и начальные приближения принимались для всех методов равными. Результаты поиска корня х= 3, представленные на рисунках, таковы. Наиболее медленно сходится метод дихотомии - 22 итерации, самый быстрый - метод простых итераций при b = -0.2 - 5 итераций. Здесь нет противоречия с утверждением, что метод Ньютона является самым быстрым.
Производная исследуемой функции в точке х = 3 равна -0.2, то есть расчет в данном случае велся практически методом Ньютона с величиной производной в точке корня уравнения. При изменении коэффициента b скорость сходимости падает и постепенно сходящийся процесс сначала зацикливается, потом становится расходящимся.
Достоинством итерационных методов является их применимость к плохо обусловленным системам и системам высоких порядков, их самоисправляемость и простота реализации на ПК. Итерационные методы для начала вычисления требуют задания какого-либо начального приближения к искомому решению.
Следует заметить, что условия и скорость сходимости итерационного процесса существенно зависят от свойств матрицы А системы и от выбора начальных приближений.
Для применения метода итераций исходную систему (2.1) или (2.2) необходимо привести к виду
после чего итерационный процесс выполняется по рекуррентным формулам
, k = 0, 1, 2, ... . (2.26а )
Матрица G и вектор получены в результате преобразования системы (2.1).
Для сходимости (2.26а ) необходимо и достаточно, чтобы |l i (G )| < 1, где l i (G ) – все собственные значения матрицы G . Сходимость будет также и в случае, если ||G || < 1, так как |l i (G )| < " ||G ||, где " – любой.
Символ || ... || означает норму матрицы. При определении ее величины чаще всего останавливаются на проверке двух условий:
||G || = или ||G || = , (2.27)
где . Сходимость гарантирована также, если исходная матрица А имеет диагональное преобладание, т. е.
. (2.28)
Если (2.27) или (2.28) выполняется, метод итерации сходится при любом начальном приближении . Чаще всего вектор берут или нулевым, или единичным, или берут сам вектор из (2.26).
Существует много подходов к преобразованию исходной системы (2.2) с матрицей А для обеспечения вида (2.26) или выполнения условий сходимости (2.27) и (2.28).
Например, (2.26) можно получить следующим образом.
Пусть А = В + С , det В ¹ 0; тогда (B + С )= Þ B = −C + Þ Þ B –1 B = − B –1 C + B –1 , откуда= − B –1 C + B –1 .
Положив –B –1 C = G , B –1 = , получим (2.26).
Из условий сходимости (2.27) и (2.28) видно, что представление А = В + С не может быть произвольным.
Если матрица А удовлетворяет условиям (2.28), то в качестве матрицы В можно выбрать нижнюю треугольную:
, a ii ¹ 0.
; Þ ; Þ ; Þ
Подбирая параметр a, можно добиться, чтобы ||G || = ||E + aA || < 1.
Если имеет место преобладание (2.28), тогда преобразование к (2.26) можно осуществить, решая каждое i -е уравнение системы (2.1) относительно x i по следующим рекуррентным формулам:
(2.28а )
Если в матрице А нет диагонального преобладания, его нужно добиться с помощью каких-либо линейных преобразований, не нарушающих их равносильности.
В качестве примера рассмотрим систему
(2.29)
Как видно, в уравнениях (1) и (2) нет диагонального преобладания, а в (3) есть, поэтому его оставляем неизменным.
Добьемся диагонального преобладания в уравнении (1). Умножим (1) на a, (2) на b, сложим оба уравнения и в полученном уравнении выберем a и b так, чтобы имело место диагональное преобладание:
(2a + 3b) х 1 + (–1,8a + 2b) х 2 +(0,4a – 1,1b)х 3 = a .
Взяв a = b = 5, получим 25х 1 + х 2 – 3,5х 3 = 5.
Для преобразования уравнения (2) с преобладанием (1) умножим на g, (2) умножим на d и из (2) вычтем (1). Получим
(3d – 2g) х 1 + (2d + 1,8g) х 2 +(–1,1d – 0,4g) х 3 = −g .
Положив d = 2, g = 3, получим 0 х 1 + 9,4 х 2 – 3,4 х 3 = −3. В результате получим систему
(2.30)
Такой прием можно применять для нахождения решения широкого класса матриц.
или
Взяв в качестве начального приближения вектор = (0,2; –0,32; 0) Т , будем решать эту систему по технологии (2.26а ):
k = 0, 1, 2, ... .
Процесс вычисления прекращается, когда два соседних приближения вектора решения совпадают по точности, т. е.
.
Технология итерационного решения вида (2.26а ) названа методомпростой итерации .
Оценка абсолютной погрешности для метода простой итерации:
где символ || ... || означает норму.
Пример 2.1 . Методом простой итерации с точностью e = 0,001 решить систему линейных уравнений:
Число шагов, дающих ответ с точностью до e = 0,001, можно определить из соотношения
£ 0,001.
Оценим сходимость по формуле (2.27). Здесь ||G || = = max{0,56; 0,61; 0,35; 0,61} = 0,61 < 1; = 2,15. Значит, сходимость обеспечена.
В качестве начального приближения возьмем вектор свободных членов, т. е. = (2,15; –0,83; 1,16; 0,44) Т . Подставим значения вектора в (2.26а ):
Продолжив вычисления, результаты занесем в таблицу:
k | х 1 | х 2 | х 3 | х 4 |
2,15 | –0,83 | 1,16 | 0,44 | |
2,9719 | –1,0775 | 1,5093 | –0,4326 | |
3,3555 | –1,0721 | 1,5075 | –0,7317 | |
3,5017 | –1,0106 | 1,5015 | –0,8111 | |
3,5511 | –0,9277 | 1,4944 | –0,8321 | |
3,5637 | –0,9563 | 1,4834 | –0,8298 | |
3,5678 | –0,9566 | 1,4890 | –0,8332 | |
3,5760 | –0,9575 | 1,4889 | –0,8356 | |
3,5709 | –0,9573 | 1,4890 | –0,8362 | |
3,5712 | –0,9571 | 1,4889 | –0,8364 | |
3,5713 | –0,9570 | 1,4890 | –0,8364 |
Сходимость в тысячных долях имеет место уже на 10-м шаге.
Ответ : х 1 » 3,571; х 2 » –0,957; х 3 » 1,489; х 4 » –0,836.
Это решение может быть получено и с помощью формул (2.28а ).
Пример 2.2 . Для иллюстрации алгоритма с помощью формул (2.28а ) рассмотрим решение системы (только две итерации):
; . (2.31)
Преобразуем систему к виду (2.26) согласно (2.28а ):
Þ (2.32)
Возьмем начальное приближение = (0; 0; 0) Т . Тогда для k = 0 очевидно, что значение = (0,5; 0,8; 1,5) Т . Подставим эти значения в (2.32), т. е. при k = 1 получим = (1,075; 1,3; 1,175) Т .
Ошибка e 2 = = max(0,575; 0,5; 0,325) = 0,575.
Блок-схема алгоритма нахождения решения СЛАУ по методу простых итераций согласно рабочим формулам (2.28а ) представлена на рис. 2.4.
Особенностью блок-схемы является наличие следующих блоков:
– блок 13 – его назначение рассмотрено ниже;
– блок 21 – вывод результатов на экран;
– блок 22 – проверка (индикатор) сходимости.
Проведем анализ предложенной схемы на примере системы (2.31) (n = 3, w = 1, e = 0,001):
= ; .
Блок 1. Вводим исходные данные A , , w, e, n : n = 3, w = 1, e = 0,001.
Цикл I . Задаем начальные значения векторов x 0 i и х i (i = 1, 2, 3).
Блок 5. Обнуляем счетчик числа итераций.
Блок 6. Обнуляем счетчик текущей погрешности.
В цикле II выполняется изменение номеров строк матрицы А и вектора .
Цикл II: i = 1: s = b 1 = 2 (блок 8).
Переходим во вложенный цикл III, блок9 – счетчик номеров столбцов матрицы А : j = 1.
Блок 10: j = i , следовательно, возвращаемся к блоку 9 и увеличиваем j на единицу: j = 2.
В блоке 10 j ¹ i (2 ¹ 1) – выполняем переход к блоку 11.
Блок 11: s = 2 – (–1) × х 0 2 = 2 – (–1) × 0 = 2, переходим к блоку 9, в котором j увеличиваем на единицу: j = 3.
В блоке 10 условие j ¹ i выполняется, поэтому переходим к блоку 11.
Блок 11: s = 2 – (–1) × х 0 3 = 2 – (–1) × 0 = 2, после чего переходим к блоку 9, в котором j увеличиваем на единицу (j = 4). Значение j больше n (n = 3) – заканчиваем цикл и переходим к блоку 12.
Блок 12: s = s / a 11 = 2 / 4 = 0,5.
Блок 13: w = 1; s = s + 0 = 0,5.
Блок 14: d = | x i – s | = | 1 – 0,5 | = 0,5.
Блок 15: x i = 0,5 (i = 1).
Блок 16. Проверяем условие d > de : 0,5 > 0, следовательно, переходим к блоку 17, в котором присваиваем de = 0,5 и выполняем возврат по ссылке «А » к следующему шагу цикла II – к блоку7, в котором i увеличиваем на единицу.
Цикл II: i = 2: s = b 2 = 4 (блок 8).
j = 1.
Посредством блока 10 j ¹ i (1 ¹ 2) – выполняем переход к блоку 11.
Блок 11: s = 4 – 1 × 0 = 4, переходим к блоку 9, в котором j увеличиваем на единицу: j = 2.
В блоке 10 условие не выполняется, поэтому переходим к блоку 9, в котором j увеличиваем на единицу: j = 3. По аналогии переходим к блоку 11.
Блок 11: s = 4 – (–2) × 0 = 4, после чего заканчиваем цикл III и переходим к блоку 12.
Блок 12: s = s / a 22 = 4 / 5 = 0,8.
Блок 13: w = 1; s = s + 0 = 0,8.
Блок 14: d = | 1 – 0,8 | = 0,2.
Блок 15: x i = 0,8 (i = 2).
Блок 16. Проверяем условие d > de : 0,2 < 0,5; следовательно, возвращаемся по ссылке «А » к следующему шагу цикла II – к блоку7.
Цикл II: i = 3: s = b 3 = 6 (блок 8).
Переходим во вложенный цикл III, блок9: j = 1.
Блок 11: s = 6 – 1 × 0 = 6, переходим к блоку 9: j = 2.
Посредством блока 10 выполняем переход к блоку 11.
Блок 11: s = 6 – 1 × 0 = 6. Заканчиваем цикл III и переходим к блоку 12.
Блок 12: s = s / a 33 = 6 / 4 = 1,5.
Блок 13: s = 1,5.
Блок 14: d = | 1 – 1,5 | = 0,5.
Блок 15: x i = 1,5 (i = 3).
Согласно блоку 16 (с учетом ссылок «А » и «С ») выходим из цикла II и переходим к блоку18.
Блок 18. Увеличиваем число итераций it = it + 1 = 0 + 1 = 1.
В блоках 19 и 20 цикла IV заменяем начальные значения х 0 i полученными значениями х i (i = 1, 2, 3).
Блок 21. Выполняем печать промежуточных значений текущей итерации, в данном случае: = (0,5; 0,8; 1,5) T , it = 1; de = 0,5.
Переходим к циклу II на блок 7 и выполняем рассмотренные вычисления с новыми начальными значениями х 0 i (i = 1, 2, 3).
После чего получим х 1 = 1,075; х 2 = 1,3; х 3 = 1,175.
Здесь , значит, метод Зейделя сходится.
По формулам (2.33)
k | х 1 | х 2 | х 3 |
0,19 | 0,97 | –0,14 | |
0,2207 | 1,0703 | –0,1915 | |
0,2354 | 1,0988 | –0,2118 | |
0,2424 | 1,1088 | –0,2196 | |
0,2454 | 1,1124 | –0,2226 | |
0,2467 | 1,1135 | –0,2237 | |
0,2472 | 1,1143 | –0,2241 | |
0,2474 | 1,1145 | –0,2243 | |
0,2475 | 1,1145 | –0,2243 |
Ответ: x 1 = 0,248; x 2 = 1,115; x 3 = –0,224.
Замечание . Если для одной и той же системы методы простой итерации и Зейделя сходятся, то метод Зейделя предпочтительнее. Однако на практике области сходимости этих методов могут быть различными, т. е. метод простой итерации сходится, а метод Зейделя расходится и наоборот. Для обоих методов, если ||G || близка к единице , скорость сходимости очень малая.
Для ускорения сходимости используется искусственный прием – так называемый метод релаксации . Суть его заключается в том, что полученное по методу итерации очередное значение x i ( k ) пересчитывается по формуле
где w принято изменять в пределах от 0 до 2 (0 < w £ 2) с каким-либо шагом (h = 0,1 или 0,2). Параметр w подбирают так, чтобы сходимость метода достигалась за минимальное число итераций.
Релаксация – постепенное ослабление какого-либо состояния тела после прекращения действия факторов, вызвавших это состояние (физ. техн.).
Пример 2.4 . Рассмотрим результат пятой итерации с применением формулы релаксации. Возьмем w = 1,5:
Как видно, получен результат почти седьмой итерации.