Итерационные методы решения нелинейных уравнений

ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ Основные этапы решения нелинейных уравнений Определение 2.1. Нелинейным уравнением называется уравнение вида f x   0 , где f  x   нелинейная функция вида: – нелинейная алгебраическая an x  an 1x n n 1 функция (полином или (2.1) многочлен)  ...  a1 x  a0 ; – трансцендентная функция – тригонометрическая, обратная тригонометрическая, логарифмическая, показательная, гиперболическая функция; –   комбинирование этих функций, например x  sin x . 2 Определение 2.2. Решением нелинейного уравнения (2.1) называется такое значение ** x , которое при подстановке в уравнение (2.1) обращает его в тождество. На практике не всегда удается найти точное решение. В этом случае решение уравнения (2.1) находят с применением приближенных (численных) методов. Определение 2.3. Приближенным решением нелинейного уравнения (2.1) называется * такое значение x , при подстановке которого в уравнение (2.1) последнее будет выполнять-     , где   малая положительная вели- ся с определенной степенью точности, т.е. f x * чина. Нахождение приближенных решений составляет основу численных методов и вычислительной математики. Решение нелинейных уравнений распадается на два этапа: отделение корней уравнений и уточнение корней нелинейных уравнений. На первом этапе необходимо исследовать уравнение и выяснить, имеются корни или нет. Если корни имеются, то необходимо определить их количество и затем найти интервалы, в каждом из которых находится только один корень, т.е. отделить корни. y y  1 x  y y  2 x   Рис. 2.1 x Первый способ отделения корней – графический. Данный метод позволяет определить количество корней на отрезке, но не единственность корня. Если f  x  имеет простой аналитический вид, то, исходя из уравнения (2.1), можно построить график функции y  f  x  . Тогда точки пересечения графика функции с осью абсцисс будут являться приближенными значениями корней исходного нелинейного уравнения. Если f  x  имеет сложный аналитический вид, то можно представить ее в виде раз- ности двух более простых функций f  x   1  x    2  x  . Так как f  x   0 , то выполняет1 ся равенство 1  x    2  x  . Построим два графика y1  1  x  , y 2   2  x  (рис. 2.1). Тогда задача решения нелинейного уравнения (2.1) сводится к поиску абсцисс точек пересечения двух графиков, которые и будут являться приближенными значениями корней уравнения (2.1). x Пример 2.1. Пусть дано нелинейное уравнение вида x  e  0 . Для решения его графическим методом представим уравнение (2.1) в виде 1  x    2  x   0 , где 1  x   x ; 2 x   e  x . Графики функций y  x ; y  e  x представлены на рис. 2.2, из которого видно, что исходное уравнение имеет единственный корень  . yx y 1 y  ex  x Рис. 2.2. Пример 2.2. Пусть задано нелинейное уравнение вида e Построив два графика функций y   x и y  e нение не имеет корней (рис. 2.3). x y  e x y  x x Рис. 2.3 2  x  0 или  x  e  x . , нетрудно заметить, что исходное урав- y 1 x Пример 2.3. Для нелинейного уравнения вида x  sin 2 x  0 с помощью аналогичных преобразований получим, что исходное уравнение имеет три корня (рис. 2.4). y yx y  sin 2 x 1   2 2 3  2  x Рис. 2.4 Второй способ отделения корней нелинейных уравнений – аналитический. Процесс отделения корней здесь основывается на следующих теоремах. Теорема 2.1. Если функция f  x  непрерывна на отрезке a, b и на концах отрезка принимает значения разных знаков (т.е. f a  f b   0 ), то на a, b содержится хотя бы один корень. Теорема 2.2. Если функция f  x  непрерывна на отрезке a, b, выполняется условие вида f a  f b   0 и производная f  x  сохраняет знак на a, b, то на отрезке имеется единственный корень. Теорема 2.3. Если функция f  x  является многочленом n -й степени и на концах от- резка a, b принимает значения разных знаков, то на a, b имеется нечетное количество корней. Если на концах отрезка a, b функция не меняет знак, то уравнение (2.1) либо не имеет корней на a, b, либо имеет четное количество корней. При аналитическом методе исследований необходимо выявить интервалы монотонности функции f  x  . Для этого необходимо найти критические точки 1 ,  2 ,...,  n , т.е. точки, в которых первая производная f  i  равна нулю или не существует. Тогда вся числовая ось разбивается на интервалы монотонности i , i 1  , на каждом из которых определяется знак производной f  xi  , где xi  i , i 1  . Затем выделяются те интервалы монотонности i , i 1  , на которых функция f x  меняет знак, т.е. выполняется неравенство f i  f i 1   0 . На каждом из этих интервалов для поиска корня используются методы уточнения корней. Наиболее распространенными методами уточнения корня на отрезке являются итерационные (приближенные) методы: метод половинного деления (метод дихотомии), метод простых итераций, метод Ньютона (метод касательных) и его модификация. 3 Метод половинного деления Для уточнения корня нелинейного уравнения (2.1) на отрезке a, b, где f a  f b   0 , а производная сохраняет знак, применим метод половинного деления. Для этого разделим отрезок a, b пополам и исследуем знак функции в полученной точке с , где c ab . Из двух отрезков a, c и c, b выбираем тот, на котором функция меняет знак. 2 Уменьшая новый отрезок в два раза, повторяем процесс и т.д. Получим последовательность отрезков a1 , b1 , a2 , b2 ,..., an , bn ,... , на концах которых выполняется неравенство f an  f bn   0 (2.2) и длины этих отрезков равны 1 (2.3) b  a  . 2n Последовательность a1 , a2 ,..., an ,... является монотонной неубывающей ограниченной последовательностью; а b1 , b2 ,..., bn ,...  монотонной невозрастающей ограниченной bn  an  последовательностью. Следовательно, эти последовательности сходятся. Перейдем к пределу при n  в левой и правой частях соотношения (2.3), получим 1 b  a   0 . n  2 n lim bn  an   lim n  Тогда lim an  lim bn   . С другой стороны, из неравенства (2.2) следует, что n  n  2 lim f bn  f an    f    0 . Последнее неравенство возможно только тогда, когда n  f   0 . Следовательно,  является корнем исходного уравнения (2.1). ПРИМЕР. Решить нелинейное уравнение 5 cos( x )  3x  0 на промежутке [0.8; 1.1] с точностью 0,0001 (но не более 3 итераций) методом половинного деления. Решение: исходим из того, что должно выполняться условие f an  f bn   0 . 2 f (0.8)  5 cos(0.82 )  3  0.8  1,6105 f (1.1)  5 cos(1.12 )  3 1.1  1,5349 f (0.8) f (1.1)  1.6105 (1.5349)  0  верно! Действия по Методу: 1. Итерация 1, i=0: a0  0.8, b0  1.1 x0  a0  b0 0.8  1.1   0.95 2 2 Для проверки точности: | f ( x0 ) || 5 cos(0.952 )  3  0.95 | 0.2482  0.0001- не верно!!! Для перехода на следующую итерацию должно быть проверено условие: f (ai ) f ( xi )  0, то ai1  ai , bi1  xi f (ai ) f ( xi )  0, то ai1  xi , bi1  bi 4 f (0.8) f (0.95)  1.6105 0.2482  0, то a1  x0  0.95, b1  b0  1.1 2. Итерация 2, i=1: a1  0.95, b1  1.1 x1  a1  b1 0.95  1.1   1.025 2 2 Для проверки точности: | f ( x1 ) || 5 cos(1.0252 )  3 1.025 || 0.5899| 0.5899  0.0001- не верно!!! Определим следующие границы: f (0.95) f (1.025)  0.2482 (0.5899)  0, то a2  a1  0.95, b2  x1  1.025 3. Итерация 3, i=2: a2  0.95, b2  1.025 x2  a2  b2 0.95  1.025   0.9875 2 2 Для проверки точности: | f ( x2 ) || 5 cos(0.98752 )  3  0.9875|| 0.1573| 0.1573  0.0001- не верно!!! Определим точность: | f ( x2 ) | 0.1573 Значит   0.2 Метод простых итераций Пусть известно, что нелинейное уравнение f  x   0 , где f  x  - непрерывная функ- ция, имеет на отрезке a, b единственный вещественный корень   a, b . Требуется найти этот корень с заданной точностью  . Применяя тождественные преобразования, приведем уравнение (2.1) к виду (2.4) x  x . Выберем произвольно приближенное значение корня (начальное приближение) x0  a, b и вычислим первое приближение  x0   x1 . Найденное значение x1 подставим в правую часть соотношения (2.4) и вычислим  x1   x2 , и так далее, т.е. xn1  xn , n  0, 1, 2,... (2.5) Продолжая процесс вычислений дальше, получим числовую последовательность x0 , x1 , x2 ,.. . Если существует предел этой последовательности, то он и является приближенным значением корня уравнения (2.4 Следовательно, предел последовательности xn  является корнем уравнения (2.4). Таким образом, корень можно вычислить с заданной точностью по следующей итера5 ционной формуле xn 1  xn , n  0,1,2,... Геометрическая интерпретация метода простых итераций Если на интервале [a,b] функция возрастает, то f ’(x)>0, если убывает - f ’(x)<0. Точками экстремума являются те точки, на которых производная меняет свой знак, при этом, если с «+» на «-», то это точки максимума, наоборот – минимума. Геометрически метод итерации может быть пояснен следующим образом. Построим на плоскости XOY графики функций y  x и y   x  . Действительный корень  уравнения (2.4) является абсциссой точки пересечения кривой y   x  с прямой y  x (рис. 2.5). Начиная процесс с некоторой точки B0 x0 , x0  , строим ломаную линию y yx y  x  A0 B0 A1 B1 B2 0 a  x2 x1 x0 b x Рис. 2.5 B0 A0 B1 A1B2 ... («лестница»), звенья которой попеременно параллельны оси OX и оси OY , вершины B0 , B1 , B2 ... лежат на кривой y  x  , а вершины A0 , A1 ,...  на прямой y  x . Общие абсциссы точек A0 и B1 , A1 и B2 , … представляют собой соответственно последовательные приближения x1 , x2 ,... корня  . В рассмотренном случае кривая y  x  пологая, x   0 и  x   1. Возможен другой вид ломаной B0 A0 B1 A1 B2 ... («спираль») (рис. 2.6). В этом случае последовательные приближения x1 , x2 ,... стремятся к корню  то с одной, то с другой стороны. В этом случае x   0 , но  x   1. y  x  y yx B1 y0 A0 y1 a M A2 A1 B2 B0 x1  x2 x0 Рис. 2.6 6 b x Однако если рассмотреть случай, где  x   1 (рис. 2.7), то процесс итераций расходится, т.е. последовательные приближения x0 , x1 , x2 ,... все дальше удаляются от корня  и в какой-то момент могут выйти за пределы отрезка a, b. Поэтому для практического применения метода простых итераций нужно определить достаточные условия сходимости итерационного процесса. Достаточное условие, при котором итерационный процесс, заданный формулой (2.5), сходится, определяет следующая теорема. y  x  y B2 yx B1 A1 B0 A0 0 a  x0 x1 b x2 x Рис. 2.7 Теорема 2.4. Пусть функция x  определена и дифференцируема на отрезке a, b, причем все ее значения x [a, b] и выполняется условие  x   q  1 при a  x  b , (2.6) тогда процесс итераций, определяемый формулой (2.5), сходится независимо от выбора начального приближения x0  a, b и предельное значение   lim xn является единственn  ным корнем уравнения (2.4) на отрезке a, b. Приведение нелинейного уравнения f  x   0 к виду x  x  , допускающему сходящиеся итерации Выполнения достаточного условия сходимости можно добиться путем перехода от исходного уравнения f  x   0 к эквивалентному виду x  x  следующим образом: умножим обе части уравнения (2.1) на неизвестную постоянную c  const  0 , c  1, затем прибавим к обеим частям переменную x , тогда получим x  cf  x   x . Обозначим через x   x  cf x  , тогда x  x . Константа c выбирается так, чтобы выполнялось достаточное условие сходимости итерационного процесса (2.6), т.е.  x   1  cf  x   1 для всех x  a, b. Это условие равносильно условию  1  1  cf  x   1, отсюда: 2  c  0 при f x   0, x  a, b ; 1) f  x  2 2) 0  c  при f  x   0, x  a, b . f  x  7 Условия окончания итерационного процесса Процесс итераций заканчивается при одновременном выполнении двух условий: 1) если два последующих приближения отличаются между собой по модулю на величину, не превышающую заданной точности  , т.е. xn 1  xn   . Отдельно этого критерия недостаточно, так как в случае крутизны графика, данное условие будет выполнено, но xn 1 может находиться далеко от корня; 2) мера удовлетворения уравнению (2.1) последнего приближения корня: f xn 1    . Отдельно данного критерия недостаточно, так как при пологой функции f  x  это условие может быть выполнено, но xn 1 может находиться далеко от корня. Метод простых итераций имеет два достоинства:  является универсальным, простым для реализации на ЭВМ и самоисправляющимся, т.е. любая неточность на каком - либо шаге итераций не отразится на конечном результате, а отразится лишь на количестве итераций. Подобные ошибки устойчивы даже по отношению к грубым ошибкам (сбоям ЭВМ), если только ошибка не выбрасывает очередное приближение за пределы области сходимости;  позволяет достигнуть любой заданной точности при любом начальном приближении x0  a, b . Недостатки метода:  трудоемкость процесса приведения уравнения (2.1) к виду (2.4);  если начальное приближение x0 выбрано достаточно далеко от корня, то число итераций, необходимых для достижения заданной точности, будет достаточно большое и объем вычислений возрастет. Пример 2.4. Доказать графическим и аналитическим методами существование единственного корня нелинейного уравнения f ( x)  e x  x  0 (2.7) на отрезке x  [1,0] и построить рабочие формулы метода простых итераций для поиска корня. 1. Докажем графическим методом единственность корня нелинейного уравнения (2.14). Из графика функции f ( x)  e  x на рис. 2.8 видно, что функция f (x) пересекает ось OX в одной точке, являющейся приближенным значением корня нелинейного уравнения (2.7). Но так как данная функция имеет сложный аналитический вид, то преобразуем x уравнение (2.7) к виду e   x и построим два графика функций y  e и y   x , имеющих более простой аналитический вид (рис. 2.9). Абсцисса точки пересечения графиков является приближенным значением корня. x x 8 y y  x 1,200 y  ex 1,000 y 0,800 0,600 0,400 0,200 0,000 -0,200 -1 -0,9 -0,8 -0,7 -0,6 -0,5 -0,4 -0,3 -0,2 -0,1 x  0 Рис. y 2.9 -0,400 -0,600 -0,800 Рис. 2.8 Для доказательства единственности корня на отрезке воспользуемся аналитическим методом. Функция f ( x) непрерывна на отрезке [ 1,0] , имеет на концах отрезка разные знаки ( f (1)  0.632; f (0)  1), а производная функции f ( x) не меняет знак на отрезке ( f ( x)  e  1  0 x [1,0] ). Следовательно, нелинейное уравнение (2.7) имеет на указанном отрезке единственный корень. x 1) Функция y=f(x) выпукла вниз на промежутках, где вторая производная положи- тельна f ' ' ( x)  0 . 3) Функция y=f(x) выпукла вверх на промежутках, где вторая производная отрицательна f ' ' ( x)  0 . 3) Функция y=f(x) имеет критические точки второго рода в точках, в которых вторая производная равна нулю или не существует (речь идет только о внутренних точках области определения функции. Точки на концах области определения не рассматриваем). 4) Функция y=f(x) имеет точки перегиба в точках, в которых вторая производная меняет знак. Метод Ньютона (метод касательных) Пусть известно, что нелинейное уравнение f  x   0 имеет на отрезке  a, b  един- ственный вещественный корень   a, b . Причем, производные f   x  , f   x  – непрерывны и сохраняют определенные знаки на отрезке  a, b  . Требуется найти этот корень с заданной точностью  . Найдем какое-либо n -е приближенное значение корня xn   ( a  xn  b ) и уточним его методом Ньютона следующим образом. Пусть   xn   n . 9 (2.8) По формуле Тейлора получим 0  f   f xn   n   f xn    n f xn  . f  xn  Следовательно,  n   . f  xn  Внося эту правку в формулу (2.8), получим рабочую формулу метода Ньютона вида: xn 1  xn  f  xn  , n  0,1, f   xn  (2.9) Геометрически метод Ньютона эквивалентен замене небольшой дуги кривой y  f  x  касательной, проведенной в некоторой точке  xn , yn  этой кривой. Для определенности положим f   x   0 и f  b   0 . Выберем начальное прибли- жение x0  b , для которого f  b   0 . Проведем касательную к кривой y  f  x  в точке B0  x0 , f  x0   . За первое приближение x1 берем точку пересечения касательной с осью OX . На кривой определим точку B1  x1, f  x1   и проведем касательную к кривой y  f  x  в этой точке. Найдем следующее приближение x2 и т.д. (рис. 2.10). y B0 y  f  x B1 B2 a  x0  x2 x1 x0  b x1 x Рис. 2.10 Если в качестве начального приближения взять другой конец отрезка  a, b  x0  a , то следующее приближение x1   a, b  . Теорема 2.5. Если f  a  f  b   0 и производные f   x  , f   x  не равны нулю и сохраняют определенные знаки на отрезке [ a, b] , то исходя из начального приближения x0 , удовлетворяющего неравенству f  x0  f   x0   0 , по методу Ньютона, заданному форму- лой (2.9), можно вычислить единственный корень  уравнения (2.1) с любой степенью точности. Замечание. Чем больше числовое значение f   x  в окрестности корня  , тем меньше правка  n . Поэтому методом Ньютона удобно пользоваться, когда в окрестности искомого корня  график функции y  f  x  имеет большую крутизну (т.е. f   xn   , тогда n 10 f  xn   0 ). Если кривая y  f  x  вблизи точки пересечения с осью OX почти горизонf   xn  n тальная (т.е. f   xn   0 , тогда n f  xn   0 ), то применять метод Ньютона для решения f   xn  n уравнения (2.1) не рекомендуется. Достоинства метода Ньютона:  обладает достаточно большой скоростью сходимости, близкой к квадратичной;  достаточно простое получение итерационной формулы (2.17). Недостатки метода Ньютона:  сходится не при любом выборе начального приближения x0 ;  применим только, когда f   x   0 для любого x   a, b  . ПРИМЕР. Решить нелинейное уравнение 5 cos( x )  3x  0 на промежутке [0.8; 1.1] с точностью 0,0001 (но не более 3 итераций) с помощью метода Ньютона. Решение: 1. Определяем тот конец отрезка, на котором выполняется условие на сходимость: f  x0  f   x0   0 . Как видно из формулы нам необходимо вычислить значение функции на концах отрезка, а также значение второй производной. 2. Определим формулу для вычисления второй производной: 2 f ' ( x)  5( sin( x 2 )  2 x  3  10  x  sin( x 2 )  3 f ' ' ( x)  10(sin( x 2 )  cos( x 2 )  2 x  x)  10 sin( x 2 )  20x 2 cos( x 2 ) 3. Определим знак условия на сходимость: f ( x0 )  f (0.8)  5 cos(0.82 )  3  0.8  1,6105 А) f ' ' ( x0 )  f ' ' (0.8)  10 sin(0.8 )  20  0.8 cos(0.8 )  16,2388 2 2 2 f (0.8)  f ' ' (0.8)  1,6105 (16,2388)  0  неверно!!! f ( x0 )  f (1.1)  5 cos(1.12 )  3 1.1  1,5349 Б) f ' ' ( x0 )  f ' ' (1.1)  10 sin(1.1 )  20  1.1 cos(1.1 )  17,8992 2 2 2 f (1.1)  f ' ' (1.1)  1,5349 (17,8992)  0  верно!!! 4. Действия по методу: Определение первого приближения начнем с правого конца интервала. 4.1. Итерация 1, i=0: Используем формулу: f  xn  , n  0,1, f   xn  f ( x0 )  1.5349 x1  x0   1.1   0.9845. f ' ( x0 )  13.2918 xn 1  xn  Проверим достигнутую точность: 11 f ( x1 )  5 cos(0.98452 )  3  0.9845  0.1239 | f ( x1 ) | 0.1239  0.0001  неверно!!! 4.2. Итерация 2, i=1: Используем формулу: x2  x1  f ( x1 )  0.1239  0.9845   0.9734. f ' ( x1 )  11.1168 Проверим достигнутую точность: f ( x2 )  5 cos(0.97342 )  3  0.9734  0.0017 | f ( x2 ) | 0.0017  0.0001  неверно!!! 4.3. Итерация 3, i=2: Используем формулу: x3  x2  f ( x2 )  0.0017  0.9734   0.9732. f ' ( x2 )  10.9036 Проверим достигнутую точность: f ( x3 )  5 cos(0.97322 )  3  0.9732  0.0005 | f ( x3 ) | 0.0005  0.0001  неверно!!! Требуемая точность корня уравнения не достигнута, но 3 итерации были выполнены. Поэтому оценим достигнутую точность для приближения x3 : | f ( x3 ) | 0.0005  0.001. Модифицированный метод Ньютона Если производная f   x  мало изменяется на отрезке [ a, b] , то можно считать, что f   xn   f   x0  . Заменив в формуле (2.9) f   xn  на f   x0  , получим рабочую формулу модифицированного метода Ньютона: xn 1  xn  f  xn  , n  0,1,... f   x0  12 (2.10) В отличие от метода Ньютона, в модифицированном методе касательная заменяется   на прямые, параллельные касательной, проведенной в точке B0 x0 , f  x0  (рис. 2.11). y f  x0  B0 y  f  x B1 B2 a  x3 x2 x1 x0  b x Рис. 2.11 Рабочая формула метода Ньютона (2.9) для данной задачи запишется так: xn 1  xn  e xn  xn , n  0,1,2,... e xn  1 Рабочая формула модифицированного метода Ньютона (2.10) для данной задачи запишется в виде: e xn  xn xn 1  xn  x , n  0,1,2,... e 0 1 13 ЛЕКЦИЯ №4 ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ. Пусть дана система линейных алгебраических уравнений (СЛАУ) вида a11x1  a12 x2  ...  a1n xn  b1, a x  a x  ...  a x  b ,  21 1 22 2 2n n 2  ... an1x1  a2n x2  ...  ann xn  bn . (4.1) Вводя в рассмотрение матрицы  a11 a12 a a A   21 22  ... ...   an1 an 2 ... a1n   x1   b1      ... a2 n  x2  b  , x , B 2  ...   ...  ... ...       ... ann   xn   bn  (4.2) систему (4.1) можно записать в виде матричного уравнения Ax  B . Для решения систем линейных (4.3) алгебраических уравнений существуют точные методы: метод Гаусса, с помощью обратной матрицы (матричный метод), по формулам Крамера. Однако, при большом числе неизвестных применение точных методов решения затруднено. В этом случае для нахождения корней системы (4.1) целесообразнее пользоваться приближенными (численными) методами. 1 Метод простых итераций для решения систем линейных алгебраических уравнений. Пусть дана система линейных алгебраических уравнений вида (4.1). Предположим, что диагональные элементы матрицы A не равны нулю, т.е. aii  0, i  1, n (в случае равенства одного или нескольких из них нулю, с помощью перестановки уравнений или других эквивалентных преобразований можно добиться, чтобы они были отличны от нуля). Разделив i -ое уравнение системы на aii , получим:  x1  1  12 x2  13 x3  ...  1n xn ,  x     x   x  ...   x ,  2 2 21 1 23 3 2n n  .............................................................  xn   n   n1x1   n 2 x2  ...   nn 1xn 1, где коэффициенты i  (4.4) aij bi , ij   при i  j,ii  0 . aii aii Введем обозначения:  1   11 12      22 2     ,    21  ...   ... ...      n1  n 2  n  ... 1n  ...  2 n  ... ...   ...  nn  (4.5) Тогда система (4.4) примет вид: x   x Систему приближений. (4.6) будем Выбираем решать начальное (4.6) методом последовательных приближение x    ; далее вычисляем следующие приближения: 1 2 1 k 1 k x      x  , x      x  ,…, x      x  , … Если последовательность приближений 2 (4.7) 1 2 k x  , x  , x  ,..., x  ,... является сходящейся, т.е. у нее существует предел   lim x k  , то этот k  предел  является решением системы (4.6). Действительно,     lim x k 1  lim    x k      . k  k  Получили      , т.е.  – является решением системы (4.6), а система (4.6) получена из системы (4.1), следовательно,  будет являться решением исходной системы (4.1). Теорема 4.1 (достаточное условие сходимости итерационного процесса). Если для приведенной системы x     x выполнено хотя бы одно из условий: n а)  ij  1, i  1, n; j 1 n б)  ij  1, j  1, n , i 1 то процесс итерации, заданный формулой x k 1 k     x  , сходится к единственному решению этой системы, независимо от выбора начального приближения. В частности процесс итерации заведомо сходится, если элементы  ij 1 приведенной системы (4.4) удовлетворяют неравенству  ij  , где n n число неизвестных системы. Следствие. Для исходной системы (4.1) итерационный процесс сходится, если n выполнены неравенства aii   aij , i  1, n (то есть модули диагональных i  j 1 3 коэффициентов для каждого уравнения системы больше суммы модулей всех остальных коэффициентов, не считая свободных членов). Теорема 4.2 (необходимое и достаточное условие сходимости процесса итерации для системы линейных уравнений). Для сходимости процесса итераций x k 1 k     x  при любом выборе начального приближения x  и любом свободном члене  необходимо и достаточно, чтобы собственные значения матрицы  (т.е. корни характеристического уравнения det    E   0 ) были по модулю меньше единицы. ПРИМЕР 4.1. Пусть известна матрица коэффициентов при неизвестных СЛАУ, а также задана единичная матрица Е: ( ) ( ). Тогда получим, ( ) ( ) ( ( ) ) ) ). Вычислим определитель данной матрицы ( ( ( ( : ( . Приравняв полученное выражение к 0 найдем собственные значения матрицы : D=0,36. . Таким образом условия теоремы 4.2 выполнены: 4 , Приведение исходной системы к виду, удобному для применения сходящегося итерационного процесса. Теорема сходимости итерационного процесса накладывает жесткие условия на коэффициенты системы (4.3). Однако, если det A  0 , то с помощью элементарных преобразований системы (4.3) ее можно заменить эквивалентной x   x , системой такой, что условия теоремы сходимости будут выполнены. Умножим левую и правую часть соотношения (4.3) слева на матрицу A 1      , где    ij , i, j  1, n - матрица с малыми по модулю элементами. Проведем преобразования: A 1      Ax  A1   B,   x   Ax  A1   B,  Если обозначить    A     B ,    A , то получим x     x . x  A1   B   Ax. 1 Если элементы матрицы  достаточно малы по модулю, т.е.  ii  0 , то элементы матрицы  будут удовлетворять достаточному условию сходимости итерационного процесса. Итерационный процесс заканчивается, если для двух приближений x k 1 n k k 1 k k 1 k и x  выполнено условие x   x    xj   xj    , где  j 1 заданная точность. 5 Метод Зейделя. Метод Зейделя представляет собой некоторую модификацию метода простых итераций. Основная его идея состоит в том, что при вычислении  k  1 -го ранее приближения неизвестной xi учитываются уже вычисленные  k  1 -е приближения неизвестных x0 , x1,..., xi 1 . Т.е. найденное Рис.4.1  k  1 -е приближение сразу же используется для получения  k  1 -го приближения последующих координат (Рис.4.1). Предполагая, что k -е приближения xi k  корней системы (4.4) известны,  k  1 -е приближения корней будут находиться по следующим итерационным формулам метода Зейделя: (4.8) 6 Теорема 4.3 (достаточное условие сходимости метода Зейделя). Если для приведенной системы x     x выполнено хотя бы одно из условий: 1)  n  1 , где  m  max  ij ; m 1i  n j 1 n 2)  l  1 , где  l  max  ij ; 1 j  n i 1 3)  k  1 , где  k  n 2  ij , i , j 1 то процесс Зейделя сходится к единственному решению системы при любом выборе начального вектора. Запишем систему (4.8) в сокращенном виде: xi k 1 i 1  i   ij xj k 1 j 1 n k   ij xj  , i  1, n (4.9) j i   Введем обозначения:   ij  B  C , где  0  B   21  ... ...    n1  n 2 ... 0 0 ... 0 0  , ... 0 0  ...  nn 1 0  . Тогда формулу (4.9) можем переписать в матричном виде: x k 1    Bx где 7 k 1 k  Cx  , (4.10)  x k    x k 1  1  1     1    k k  1     x  x      2  , x k    2  , x k 1   2  .  ...   ...   ...        x k    x k 1   n   n   n  Теорема 4.4 (необходимые и достаточные условия сходимости метода Зейделя). Для сходимости процесса Зейделя, заданного формулой (4.9), для приведенной системы линейных уравнений (4.4) при любом выборе свободного члена  и начального вектора x  необходимо и достаточно, чтобы все корни 1, 2 ,..., n уравнения det  C   E  B     0 были по модулю меньше единицы. Пример 4.2. Построить рабочие формулы метода простых итераций и метода Зейделя для численного решения СЛАУ вида: 8 x  5 y  z  1;   x  6 y  2 z  7;  x  y  4 z  9.  (4.11) Решение. Заметим, что система (4.11) имеет точное решение x  1; y  2; z  3 . Из системы (4.11) видно, что модули диагональных коэффициентов в каждом уравнении отличны от нуля и больше суммы модулей всех остальных коэффициентов, не считая столбца свободных членов. Тогда, разделим каждое уравнение системы (4.11) на соответствующий диагональный коэффициент, сформируем столбец ( x, y, z) в левой части и перенесем остальные слагаемые в правую часть, получим рабочие формулы метода простых итераций вида: 8   k 1 1 5  k  1  k    y  z ; x 8 8 8    k 1 7 1  k  1  k    x  z ; y 6 6 3    k 1 9 1  k  1  k    x  y , k  0,1,2,... z 4 4 4  Начальное приближение обычно выбирают равным  столбцу  1 7 9 свободных членов преобразованной системы x(0) , y (0) , z (0)   , ,  . 8 6 4 Рабочие формулы метода Зейделя запишутся так:   k 1 1 5  k  1  k    y  z ; x 8 8 8    k 1 7 1  k 1 1  k    x  z ; y 6 6 3    k 1 9 1  k 1 1  k 1   x  y , k  0,1,2,... z 4 4 4  ЗАДАНИЕ. Построить рабочие формулы метода простых итераций и метода Зейделя для численного решения СЛАУ вида: { . Метод релаксации. Рассмотрим систему линейных алгебраических уравнений (4.1) a11x1  a12 x2  ...  a1n xn  b1, a x  a x  ...  a x  b ,  21 1 22 2 2n n 2  ... an1x1  a2n x2  ...  ann xn  bn , в которой aii  0, i  1, n . Сделаем преобразования: свободные члены перенесем в левую часть 9 и каждое i -ое уравнение поделим на  1 aii , i  1, n . Таким образом, получим систему, удобную для релаксации:  x1  b12 x2  ...  b1n xn  c1  0, b x  x  ...  b x  c  0,  21 1 2 2n n 2  ... bn1x1  bn 2 x2  ...  xn  cn  0. где bij  aij aii , i  j , ci  (4.12) bi . aii Введем понятие невязки для приближенного решения x0 . Пусть дана система Ax  B , тогда точное решение x можно записать  1    в виде x  x0   , где    ...  -правка корня x0 . Подставим x  x0   в    n систему, получим A  x0     B, Ax0  A  B, A  B  Ax0 . Введем обозначение B  Ax0   . Тогда A   . Выражение   B  Ax0 называется невязкой для приближенного решения x0 . Пусть задано начальное приближение системы (4.12):   x   x1  , x2  ,..., xn  . Подставим данное приближение в систему (4.12) и получим невязки Ri  , i  1, n : 10 n R1   c1  x1    b1 j xj  , j 2 R2   c2  x2   0  b2 j x j , n j 1, j  2 (4.13) ... n 1 Rn   cn  xn    bnj xj  . j 1 Если одной из неизвестных xk  дать приращение xk  , то соответствующая невязка Rk  уменьшится на величину xk  , а все остальные невязки Rq  , q  k изменятся на величину bqk xk  . Чтобы 1 обратить очередную невязку Rk  в нуль, нужно величине xk  дать 1 приращение xk   Rk  , следовательно, Rk   0 , а остальные невязки будут равны 1 Rq   Rq   bqk xk  , q  k . Метод релаксации (метод ослабления) заключается в том, что на каждом шаге обращают в нуль максимальную по модулю невязку путем изменения значения соответствующей компоненты приближения. Процесс заканчивается, когда все невязки последней преобразованной системы будут равны нулю с заданной точностью. Пример 4.3. Решить систему методом релаксации, производя вычисления с двумя десятичными знаками. 10 x1  2 x2  2 x3  6,   x1  10 x2  2 x3  7,  x  x  10 x  8. 3  1 2 Решение. Приведем систему (4.14) к виду, удобному для релаксации 11 (4.14)  x1  0.2 x2  0.2 x3  0.6  0,  0.1x1  x2  0.2 x3  0.7  0, 0.1x  0.1x  x  0.8  0. 1 2 3  (4.15) В качестве начального приближения выбираем x1   x2   x3   0 . Находим соответствующие невязки: R1   0.60 , R2   0.70 , R3   0.80 . Выбираем максимальную невязку и полагаем x3   0.80 , тогда R3   0, 1 1 R1   R1   b13x3   0.6  0.2  0.80  0.76, 1 R2   R2   b23x3   0.7  0.2  0.80  0.86. Опять выбираем максимальную невязку и полагаем x2   0.86 , тогда 1 2 R2   0, 2 2 2 R1   R1   b12x2   0.76  0.2  0.86  0.93, 2 1 1 R3   R3   b32x2   0  0.1  0.86  0.09. 2 Далее x1   0.93 и 3 R1   0, 3 2 2 R2   R2   b21x1   0  0.1  0.93  0.09, 3 2 2 R3   R3   b31x1   0.09  0.1  0.93  0.18. 3 x3   0.18 , 4 R3   0, 4 3 3 R1   R1   b13x3   0  0.2  0.18  0.04, 4 3 3 R2   R2   b23x3   0.09  0.2  0.18  0.13. 4 x2   0.13 , 12 5 R2   0, 5 4 4 R1   R1   b12x2   0.04  0.2  0.13  0.07, 5 4 4 R3   R3   b32x2   0  0.1  0.13  0.01. 5 x1   0.07 , 6 R1   0, 6 R2   0  0.1  0.07  0.01, 6 R3   0.1  0.1  0.07  0.02. 6 x3   0.02 , 7 R3   0, 7 R1   0  0.2  0.02  0.004  0.00, 7 R2   0.01  0.2  0.02  0.01. 7 x2   0.01 , 8 R2   0, 8 R1   0  0.2  0.01  0.00, 8 R3   0  0.1  0.01  0.00. Окончательно получим: 2 5 x1  x1   x1   x1   0  0.93  0.07  1.00, 1 4 7 x2  x2   x2   x2   x2   0  0.86  0.13  0.01  1.00, 3 6 x3  x3   x3   x3   x3   0  0.80  0.18  0.02  1.00. 13 ПРИБЛИЖЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И СИСТЕМ. Дифференциальные уравнения являются основным математическим инструментом моделирования и анализа разнообразных явлений и процессов в науке и технике. Методы их решения подразделяются на два класса: 1) аналитические методы, в которых решение получается в виде аналитических функций; 2) численные (приближенные) методы, где искомые интегральные кривые получают в виде таблиц их числовых значений. Применение аналитических методов позволяет исследовать полученные решения методами математического анализа и сделать соответствующие выводы о свойствах моделируемого явления или процесса. К сожалению, с помощью таких методов можно решать достаточно ограниченный круг реальных задач. Численные методы позволяют получить с определенной точностью приближенное решение практически любой задачи. Решить дифференциальное уравнение y  f  x,y  (8.1) численным методом означает, что для заданной последовательности аргументов x0 , x1 , x2 , . . . , xn и числа y0  y  x0  , не определяя аналитического вида функции y  F  x  , найти значения y1 , y2 , . . . ,yn , удовлетворяющие условиям: F  x0   y0 , yk  F  xk  , k  1 , n . Рассмотрим три наиболее распространенных при решении практических задач численных метода интегрирования: метод Эйлера, метод РунгеКутта и метод Адамса. МЕТОД ЭЙЛЕРА Этот метод обладает малой точностью и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других численных методов. Пусть дано дифференциальное уравнение с начальным условием (задача Коши) y  f  x,y  , y0  y  x0  (8.2) и выполняются условия существования и единственности решения. Требуется найти решение y  x  задачи Коши (8.2). Выбрав шаг h - достаточно малый, равный h   b  a  n , строим систему равноотстоящих точек x0 , x1 , . . . , xn , xi  x0  ih, i  0,n . 1 Искомую интегральную кривую y  y  x  , проходящую через точку M 0  x0 , y0  , приближенно заменим ломаной Эйлера с вершинами M i  xi yi  , i  0,1,2... (Рис.8.1). y y2 M 1 y1 M0  y0 h x0 tg  y0  y1  y0 h h x1 x2 x Рис.8.1 Звено ломаной M i M i 1 , заключенное между xi и xi 1 , наклонено к оси OX под углом  . Тангенс этого угла вычисляется по формуле: y y tg i  i 1 i  y  xi   f  xi , yi  . h Сделав преобразование, получим формулу Эйлера: yi 1  yi  h f  xi , yi  , i =0, n . (8.3) Вычисление значений y1, y2 , . . . , yn осуществляется с использованием формулы (8.3) следующим образом. По заданным начальным условиям x0  a и y0 полагая i  0 в выражении (8.3) вычисляется значение y1  y0  h f  x0 , y0  . (8.4) Далее определяя значение аргумента x по формуле x1  x0  h , используя найденное значение y1 и полагая в формуле (8.3) i  1 вычисляем следующее приближенное значение интегральной кривой y  F  x  : y2  y1  h f ( x1, y1) . (8.5) Поступая аналогичным образом при i  2, n  1 определяем все yi , в том числе последнее значение остальные значения yn  yn 1  h f  xn 1, yn 1  , которое соответствует значению аргумента xn  b . Таким образом, соединяя на координатной плоскости точки  x0 , y0  ,  x1, y1  , . . . ,  xn ,yn  отрезками прямых, получаем ломаную линию с вершинами в точках M 0  x0 , y0  , M1  x1, y1  , . . . ,M n  xn , yn  . Запишем разложение yi 1 в ряд Тейлора: Как правило, шаг h выбирают таким образом, чтобы h2   , где  2 заданная точность. Метод Эйлера может быть применен к решению систем дифференциальных уравнений. Пусть задана система двух дифференциальных уравнений первого порядка  y  f1 ( x, y, z ), (8.6)   z  f 2 ( x, y, z ), с начальными условиями y  x0   y0 , z  x0   z0 . Необходимо найти решение этой задачи Коши. Проводя аналогичные рассуждения, получаем расчетные формулы вида: yi 1  yi  hf1  xi , yi , zi  , zi 1  zi  hf 2  xi , yi , zi  , (8.7) xi 1  xi  h , i  0,1,2,..., где h  шаг интегрирования. В результате применения расчетной схемы (8.7) получается приближенное представление интегральных кривых y  F1  x  и z  F2  x  в форме двух ломаных Эйлера, построенных по полученным точкам  xi , yi  ,  xi , zi  , i  0,1,2,... . Достоинством метода Эйлера является его простота и высокая скорость поиска решения. Недостатками метода Эйлера являются малая точность и систематическое накопление ошибок, так как при вычислении значений на каждом последующем шаге исходные данные не являются точными и содержат погрешности, зависящие от неточности предшествующих вычислений. МЕТОД РУНГЕ-КУТТА(Ы) Данный метод является одним из наиболее распространенных численных методов интегрирования обыкновенных дифференциальных уравнений. По сравнению с описанным выше методом Эйлера метод РунгеКутта имеет более высокую точность, но невысокую скорость поиска решения, так как метод относится к классу многошаговых методов. Пусть дано дифференциальное уравнение первого порядка y  f  x,y  , (8.8) с начальным условием y0  y  x0  . (8.9) Выберем шаг h и для краткости введем обозначения xi  x0  ih , yi  y  xi  , где i  0,1,... . Рассмотрим числа: 3 k1i  hf  xi , yi  , k  h  k2i  hf  xi  , yi  1i  , 2 2   k  h  k3i  hf  xi  , yi  2i  , 2 2   (8.10) k4i  hf  xi  h, yi  k3i  . По методу Рунге-Кутта последовательные значения yi искомой функции y определяются по формуле: 1 yi 1  yi   k1i  2k2i  2k3i  k4i  . (8.11) 6 Погрешность метода Рунге-Кутта, заданного формулой (8.11), на каждом шаге есть величина порядка h 5 (в предположении, что f  x, y   C   ). Формулу (8.11) еще называют формулой Рунге-Кутта четвертого порядка точности. Помимо формулы (8.11) существуют еще другие формулы типа Рунге-Кутта(ы) с иными порядками точности. В частности формула yi 1  yi  k2i  формула Рунге-Кутта второго порядка точности. Эта фор5 мула на каждом шаге дает погрешность порядка h 3 . Для определения правильности выбора шага h на практике обычно на каждом этапе из двух шагов применяют двойной пересчет. Исходя из текущего верного значения y  xi  , вычисляют y  xi  2h  двумя способами: вначале с шагом h , а затем с шагом 2h . Если расхождение полученных результатов не превышает допустимой погрешности, то шаг h для данного этапа выбран правильно и полученное с его помощью значение можно принять за y  xi  2h  . В противном случае шаг уменьшается в два раза. Эту вычислительную схему легко запрограммировать на ЭВМ. Метод Рунге-Кутта(ы) может быть использован и при решении систем дифференциальных уравнений. Рассмотрим задачу Коши для системы двух дифференциальных уравнений:  y  f1 ( x, y, z ),   z  f 2 ( x, y, z ), с начальными условиями y  x0   y0 , z  x0   z0 . Формулы метода Рунге-Кутта для данной системы примут вид: 4 1  k1i  2k2i  2k3i  k4i  , 6 1 zi 1  zi   m1i  2m2i  2m3i  m4i  , i  0, n. 6 yi 1  yi  где k1i  hf1  xi , yi , zi  , m1i  hf 2  xi , yi , zi  , k2i  hf1  xi  0.5h ,yi  0.5k1i ,zi  0.5m1i  , m2i  hf 2  xi  0.5h ,yi  0.5k1i ,zi  0.5m1i  , k3i  hf1  xi  0.5h ,yi  0.5k2i ,zi  0.5m2i  , m3i  hf 2  xi  0.5h ,yi  0.5k2i ,zi  0.5m2i  , k4i  hf1  xi  h ,yi  k3i ,zi  m3i  , m4i  hf 2  xi  h ,yi  k3i ,zi  m3i  . Метод Рунге-Кутта(ы) обладает значительной точностью и, несмотря на свою трудоемкость, широко используется при численном решении дифференциальных уравнений и систем. Важным преимуществом этого метода является возможность применения переменного шага, что позволяет учитывать локальные особенности искомой функции. МЕТОД АДАМСА Этот метод численного интегрирования разработан Адамсом в 1855 году по просьбе известного английского артиллериста Башфорта, который занимался баллистикой. В последствии этот метод был забыт и вновь открыт в начале XX века норвежским математиком Штермером. Популяризация метода Адамса и дальнейшее его усовершенствование связано с именем А.Н. Крылова. Пусть дано дифференциальное уравнение первого порядка y  f  x,y  , (8.12) с начальным условием y0  y  x0  . (8.13) Пусть xi , i  0,1,...  система равноотстоящих значений с шагом h и yi  y  xi  . Очевидно, что xi 1 yi   y  x  dx . (8.14) xi Запишем вторую интерполяционную формулу Ньютона с точностью 5 до разностей четвертого порядка: q  q  1 2 q  q  1 q  2  3 y  x   yi  qyi 1   yi  2   yi 3 , (8.15) 2! 3! x  xn где q  . h В формуле (8.15) функцию y заменим на производную y , получим: q  q  1 2 q  q  1 q  2  3 y  x   yi  qyi1   yi 2   yi3 . (8.16) 2! 3! Так как dx  h  dq , то подставив (8.16) в (8.14), получим: 1  q2  q 2 q3  3q 2  2q 3 yi  h   yi  qyi1   yi 2   yi3  dq . 2 6 0  После преобразований будем иметь: 1 5 3 yi 1  yi  hyi    hyi1    2  hyi 2    3  hyi3  . (8.17) 2 12 8 Формула (8.17) называется экстраполяционной формулой Адамса. Для начала итерационного процесса нужно знать начальные значения y0 , y1, y2 , y3 , так называемый начальный отрезок, который определяют, исходя из начального условия (8.13), каким-либо численным методом (например, методом Рунге-Кутта). Зная значения y0 , y1, y2 , y3 из (8.12) находят y0 , y1 , y2 , y3 и составляют таблицу разностей:   hy0  ,   hy1  ,   hy2  ,  2  hy0  ,  2  hy1  , 3  hy0  . (8.18) Дальнейшие значения yi , i  4,5,... искомого решения можно шаг за шагом вычислять по формуле Адамса (8.17), пополняя по мере необходимости таблицу разностей (8.18). Для работы на ЭВМ формулу Адамса применяют в раскрытом виде. Так как yi1  yi  yi1,  2 yi 2  yi  2 yi1  yi 2 , 3 yi3  yi  3 yi1  3 yi 2  yi3 , то после приведения подобных членов имеем: h yi 1  yi   55 yi  59 yi1  37 yi 2  9 yi3  , (8.19) 24 xi 1  xi  h. На практике шаг h выбирают так, чтобы можно было пренебречь величиной 1 3   hyi 2  . 24 6 Метод Адамса легко распространяется на системы дифференциальных уравнений. Погрешность метода Адамса имеет тот же порядок, что и метод Рунге-Кутта. 7