Численное решение нелинейных алгебраических уравнений и систем

Напомним критерий Коши сходимости числовой последовательности: последовательность {uk}, k = 0, 1, ... является сходящейся, если для любого положительного числа $$\varepsilon$$ существует номер n такой, что при всех k > N и любых натуральных p расстояние между членами последовательности uk и uk + 1 меньше $$\varepsilon,$$ т.е. $$|u_{k} - u_{k + p}\} | < \varepsilon .$$

Напомним критерий Коши для последовательности элементов метрического пространства: последовательность {uk}, k = 0, 1, ... является сходящейся, если для любого $$\varepsilon > 0$$ существует номер n такой, что при всех k > N и любом натуральном p расстояние $$\rho (u_{k}, u_{k + p}) < \varepsilon.$$

Продолжим доказательство. Переходя в последнем неравенстве к пределу при $$p \to \infty,$$ получим

Покажем, что $$\mathbf{U}$$ есть корень уравнения (5.2)

Поскольку k выбрано произвольно, а левая часть от k не зависит, то $$\rho \left[{{\mathbf{U}}, {\mathbf{F}}({\mathbf{U}})}\right] = 0,$$ или $$\mathbf{U} = \mathbf{F}(\mathbf{U}).$$

В случае скалярного уравнения имеем $$\theta = u_{k} + t(u_{k + 1} - u_{k}),$$

откуда следует условие сходимости итерационного процесса $$\max \left|{F^{\prime}(u)}\right| \le q < 1.$$ Отрезок $$\Delta$$ включает в себя всю последовательность { uk }, $$u_k \in \Delta, k = 0, 1, 2, \ldots$$

В случае решения системы нелинейных уравнений достаточным условием сходимости итерационного процесса будет $$\left\|{{\mathbf{F^{\prime}}}({\mathbf{u}})}\right\| < 1,$$ где $${\mathbf{F^{\prime}}}({\mathbf{u}})$$ — матрица Якоби.

Теорема (без доказательства.) Пусть область $$G \in L^n$$ выпуклая, $${\mathbf{u}} \in G$$ а компоненты $$F_i(\mathbf{u})$$ вектор - функции $${\mathbf{F}}({\mathbf{u}}) = (F_1, \ldots , F_n)^T$$ имеют равномерно непрерывные производные первого порядка. Положим, что норма матрицы Якоби

не превосходит некоторого числа $$0 \le q < 1,$$ т.е. $$\left\|{\mathbf{Y}}\right\| \le q < 1$$ для всех $${\mathbf{u}} \in G.$$

В этом случае отображение $${\mathbf{v}} = \mathbf{F}(\mathbf{u})$$ является сжимающим в области G, т.е. $$\rho ({\mathbf{F}}({\mathbf{u}}_1 ), {\mathbf{F}}({\mathbf{u}}_2 )) \le q \rho ({\mathbf{u}}_1, {\mathbf{u}}_2 ),$$ или $$\left\|{{\mathbf{F}}({\mathbf{u}}_1 ) - {\mathbf{F}}({\mathbf{u}}_2 )}\right\| \le q\left\|{{\mathbf{u}}_1 - {\mathbf{u}}_2} \right\|.$$

Геометрическая интерпретация метода простой итерации для скалярного случая uk + 1 = F(uk) приведена на рис. 5.1. Алгоритм метода простых итераций таков.

Метод релаксации. Без ограничения общности рассмотрим скалярный случай. Положим $$F(u) = u + \tau f(u)$$ и построим итерационный процесс $$u_{k + 1} = u_{k} + \tau f(u_{k}), u_{0} = a.$$

Тогда $$F'(u) = 1 + \tau f'(u)$$ и $$\tau$$ выбирается из условия |F'(u)| < 1, причем, чем меньше значение |F'(u)|, тем быстрее будет сходиться итерационный процесс. В частности, если положить F'(u) = 0, то, $$\tau = - [f'(u)]^{ - 1},$$ а формулы итерационного процесса будут uk + 1 = uk - [f'(uk)] - 1 f(uk), u0 = a.

5.2. Метод Ньютона

Как и выше, необходимо найти решение уравнения f(u) = 0. Пусть uk есть k приближение решения ( k итерация ). Следующее приближение ищем в виде $$u_{k + 1} = u_{k} + \Delta u_{k},$$ разложив функцию f(u) в ряд Тейлора с точностью до членов первого порядка: $$f(u_k + \Delta u_k) = f(u_k) + f^{\prime}_u (u_k) \cdot \Delta u_k + O(\Delta ^2 u_k).$$

Пренебрегая членами $$O(\Delta ^{2}u_{k}),$$ получим линеаризованное уравнение для определения $$\Delta u_{k}$$:

Это уже знакомая формула, полученная в результате оптимизации релаксационного варианта метода простой итерации.

Для системы уравнений матрица Якоби $${\mathbf{f^{\prime}}}_{\mathbf{u}} ({\mathbf{u}}_k)$$ будет

а метод Ньютона выглядит следующим образом:

где fi = fi (u1, ..., un). Приходим к СЛАУ вида $${\mathbf{A}}(\Delta {\mathbf{u}}) = - {\mathbf{f}},$$ где $$\Delta {\mathbf{u}} = {(\Delta u_1, \Delta u_2, \ldots , \Delta u_n)}^T ,$$ $${\mathbf{f}} = {(f_1, f_2, \ldots , f_n)}^T ,$$ $$\Delta {\mathbf{u}} \in L^{n}, {\mathbf{f}} \in L^{n} ,$$ $$\Delta u_i = u_i^{k + 1} - u_i^{k} ,$$ i = 1, ..., n, k = 0, 1, ....

Геометрический смысл метода Ньютона в одномерном случае проиллюстрирован на рис. 5.2. Заменим f(u) в точках uk — каждом приближении к корню — касательными. За следующее приближение по методу Ньютона примем значение u точки пересечения касательной с осью абсцисс. Метод Ньютона называют также методом линеаризации или методом касательных.

Теорема о квадратичной сходимости метода Ньютона [5.2], [5.4]

Сформулируем и докажем теорему для одномерного (скалярного) случая. Аналогичная теорема будет справедлива и для систем нелинейных уравнений.

Теорема. Пусть существуют первые две ограниченные производные f(u) и, кроме того, существует [f'u(u)] - 1 ; причем имеют место оценки $$\left|{f^{\prime\prime}{uu}}\right| \le C_2, \left|{\left[{f^{\prime}(u)}\right]^{- 1}}\right| \le C_1$$ (отображение f(u) равномерно невырождено), а начальное приближение выбирается из условия

Тогда метод Ньютона сходится с квадратичной скоростью сходимости.

Доказательство.

Разложим f(uk + 1) в ряд Тейлора в окрестности f(uk), ограничившись квадратичными членами разложения:

Здесь введено обозначение $$\Delta u_{k} = u_{k + 1} - u_{k}.$$ Переходя к абсолютной величине и учитывая, что для метода Ньютона $$f(u_{k}) + f'_{u}(u_{k}) \Delta u_{k} = 0,$$ или $$\Delta u = - [f'_{u}(u)]^{ - 1}f(u_{k}),$$ получим

так как | [f'u(u)]- 1 | < C1.

Введем в рассмотрение невязку как rk = | f(uk) | , получим $$r_{k + 1} \le Cr_k^2 ,$$ где $$C = C_2 C_1^2 .$$

Можно выписать цепочку соотношений $$r_1 \le Cr_0^2, r_2 \le Cr_1^2 \le C^3 r_0^4, r_3 \le Cr_2^2 \le C^7 r_0^8 ,$$ и т.д., в результате для невязки на k итерации получается выражение $$r_k \le C^{- 1} (Cr_0)^{2^k}.$$

Неравенства $$r_{k + 1} \le Cr_k^2$$ и $$r_k \le C^{- 1} (Cr_0)^{2^k}$$ являются определением квадратичной скорости сходимости.

Для сходимости итерационного процесса Ньютона достаточно, чтобы было выполнено условие, следующее из последнего неравенства: $$Cr_0 = C \left|f(u_0)\right| \le q < 1.$$ Отсюда следуют ограничения на начальное приближение, в частности, $$| {f(u_0 )} | \le 1/{C}.$$ Теорема доказана.

Замечание. Несложно показать, что погрешность, определяемая, как $$\varepsilon_k = \rho ({\mathbf{u}}_k - {\mathbf{U}}),$$ или, в скалярном случае, $$\varepsilon _{k} = | u_{k} - U |,$$ убывает квадратично. Для этого разложим f(uk + 1) в окрестности uk в ряд Тейлора до первого члена (или линеаризуем f(uk + 1))

Так как в методе Ньютона приближения находятся достаточно близко к корню уравнения и $$u_{k + 1} \approx U,$$ получим

Разделив полученное равенство на f'u(uk), приходим к оценке $$$ U - \left[{u_k - \frac{f(u_k)}{f^{\prime}(u_k)}}\right] = \frac{{f^{\prime\prime}_u (\theta )}}{{2f^{\prime}_u (u_k)}}(U - u_k)^2 $,$$ откуда следует $$$ \left| U - (u_k - \frac{f (u_k)}{{f^{\prime}}_u (u_k)})\right| \le \left| \frac{\max \left|f^{\prime\prime}_u(\Theta)\right|}{2f^{\prime}_u(u_k)}\right| \left|U - u_k\right| ^2, u_k \in \Delta , k = 0, 1, \ldots $.$$ Левая часть последнего неравенства по формуле Ньютона равна $$\varepsilon _{k + 1} = | U - u_{k + 1} |.$$ В таком случае $$\varepsilon_{{k + 1}} \le \bar {C}\varepsilon _k^2,$$ где $$\bar{C} = \frac{{\max\limits_\Delta } {\left|{f^{\prime\prime}(\theta )}\right|}}{2\left|{f^{\prime}(u_k)}\right|},$$ откуда последовательно находим q(k)[f(u)].

Итерационные процессы, имеющие третий и четвертый порядок сходимости, представляются формулами

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

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

Этот метод оказывается приемлемым, поскольку начальное приближение в методе Ньютона обычно выбирается достаточно близким к корню уравнения.

Методом секущих (или разностным методом Ньютона ) называется итерационный метод, в котором вместо производной вычисляется разностное выражение $$$ f^{\prime}(u_k) \approx \frac{{f(u_k) - f(u_{k - 1})}}{{\tau_k}} $,$$ откуда $$$ u_{{k + 1}} = u_k - \frac{{f(u_k)\tau_k}}{{f(u_k) - f(u_{k - 1})}}, \tau_k = u_k - u_{k - 1} $.$$

5.3. О вариационных подходах к решению нелинейных систем уравнений

Рассмотрим систему нелинейных уравнений f(u, v) = 0, g(u, v) = 0, $${u, v} \in \Omega.$$

Рассмотрим функционал $$\Phi (u, v) = f^{2}(u, v) + g^{2}(u, v).$$

Так как $$\Phi$$ неотрицателен, то найдется точка $$\left\{{\bar {u}, \bar {v}}\right\}$$ такая, что $$\left\{{\bar {u}, \bar {v}}\right\} = \mathop {\arg }\limits_{u, v \in \Omega } \min \Phi (u, v),$$ но $$min \Phi (u, v),$$ очевидно, достигается при f(u, v) = 0, g(u, v) = 0, т.е. на решении исходной системы уравнений.

Построим итерационный процесс, соответствующий методу градиентного спуска

где $$\tau _{k}$$ — параметр, который выбирается, например, из условия минимальности $$\Phi (u_{k + 1}, v_{k + 1})$$ в данном направлении (метод наискорейшего спуска), {pk, qk} — вектор, определяющий направление минимизации. На каждом шаге итераций решается задача минимизации $$\Phi$$ по одному аргументу.

5.4. Метод Чебышёва построения итерационных процессов высшего порядка

Предположим, что существует функция g(u), обратная к f(u). При этом u = g[f(u)], U = g(0). Пусть, кроме того, f(u) непрерывна и имеет необходимое число непрерывных производных на отрезке, внутри которого лежат все члены последовательности {uk}, k = 0, 1, ... Обратная функция имеет такое же количество непрерывных производных, как и f(u). Разложим функцию g(f[v]) = g(h) в ряд Тейлора в окрестности корня - точки w = f(u)

Тогда, учитывая, что u = g[f(u)], w = f(u), h = f(v), получим

Можно показать, что итерационный метод

имеет порядок сходимости n + 1. Для вычисления производных обратной функции u = g[f(u)] воспользуемся правилом дифференцирования сложной функции:

5.5. Разностные отображения в нелинейной динамике

Рассмотрим последовательность чисел $$u_{k + 1} \in R$$ ( R — множество вещественных чисел), каждый член которой связан с предыдущим рекуррентным соотношением

где $$k \in N$$ ( N — множество натуральных чисел). Соотношения (5.4) называются разностными отображениями (уравнениями) с дискретным аргументом.

Такие уравнения появляются при моделировании процессов, в которых величина u рассматривается через определенные промежутки времени. Например, еще в середине XIX века Ферхюльст для описания динамики популяционной системы предложил измерять ежегодно численность особей uk, где k — номер года. Относительная численность uk + 1 полагалась пропорциональной численности в k год, однако она начинает убывать, когда животных становится много ( uk сравнимо с 1):

где

Другой пример из экономической области — задача о банковских сбережениях. Пусть u0 — денежный вклад, растущий в соответствии с постоянным процентом $$\delta,$$ по закону:

Пусть далее законодательный орган, желая воспрепятствовать такому обогащению вкладчика, издает закон о том, чтобы процент убывал пропорционально uk, т.е.

Тогда счет в банке изменился бы по закону

т.е. в соответствии с моделью (5.6).

Так как $$u_k \in [0, 1],$$ то $$\lambda \in [0, 4].$$ Отображение (5.6) называется логистическим. К нему можно также придти, применив простейший из численных методов решения обыкновенных дифференциальных уравнений — явный метод Эйлера (лекция 8) для решения дифференциального уравнения динамики популяции (уравнения Ферхюльста)

где u — численность популяции. Вводя шаг по времени, получим разностный аналог уравнения (5.8):

откуда получаем

И отображение (5.7), и отображение (5.10) легко приводится к виду (5.6). Достаточно произвести очевидную замену переменных. Для (5.10) эта замена будет $$$ \mu = \alpha , z_k = \frac{\alpha }{\beta }u_k $.$$ Как двумерное обобщение логистического отображения можно рассматривать отображение Хенона

или

К довольно известным двумерным дискретным моделям относится также отображение Чирикова, предложенное для моделирования поведения незатухающего ротатора, возбуждаемого внешними толчками:

Рассмотрим подробнее свойства отображения:

Заметим, что f(0) = f(1) = 0 и $$max f(u) = f(0, 5) = \lambda / 4,$$ то при $$0 < \lambda < 4$$ интервал X = [0, 1] отображается в себя, $$u \in X.$$

Введем обозначения

Последовательность f, f2, ..., fk, ... называется траекторией отображения и обозначается $$\left\{f^{k}(u_0) \right\}_0^\infty.$$

Определение. Точка $$a \in X$$ ( X — множество, включающее в себя все значения отображения (5.4)) называется предельной точкой траектории $$\left\{f^{k} (u_0 ) \right\}_{k = 0}^\infty,$$ если существует последовательность $$k_1 < k_2 < ... < k_n \mathop \to\limits_{n \to \infty } \infty$$ такая, что $$f^{k_n} \to a, n = 1, 2, \ldots.$$

Рассмотрим вначале случай $$0 < \lambda < 1.$$ На X = [0, 1] существует только одна предельная (или неподвижная) точка x = 0. Любая последовательность, $$\left\{f^{k} (u_0 ) \right \}_{k = 0}^\infty$$ сходится к предельной точке рассматриваемого отображения x = 0. Если рассматривается популяционная модель, то это означает, что рассматриваемая популяция не может выжить.

Из теоремы о сжимающем отображении следует, что последовательность $$\left\{{u_k}\right\}_{n = 0}^\infty$$ сходится к своей предельной точке, если $$\left|{f^{\prime}_u}\right| \le 1.$$

В этом случае точка называется притягивающей. При выполнении условия | f'u | > 1 точка называется отталкивающей.

Графическое изображение траектории (лесенка Ламерея) представлено на рис. 5.3.

Теперь рассмотрим случай $$1 < \lambda < 3.$$

В случае, когда $$\lambda > 1$$ — неподвижная точка, u = 0 становится отталкивающей, поскольку | f'(0) | > 1, а на отрезке [0, 1] появляется другая неподвижная точка $$u_{1} = 1 - \lambda ^{ - 1}.$$

Производная для рассматриваемого отображения $$| f'(u_{1}) | = | 2 - \lambda | < 1.$$ Точка u1 при $$1 < \lambda \le 3$$ является притягивающей.

Отметим, что при $$1 < \lambda \le 2$$ производная f'(u1) > 0 и траектория $$\left\{f^{k} (u_0 ) \right\}_{k = 1}^\infty $$ стремится монотонно к u1 ( рис. 5.4); при $$2 < \lambda \le 3$$ производная f'(u1) < 0 и траектория приближается к u1 немонотонно, поочередно принимая значения то меньше, то больше этого значения.

При $$\lambda = 3$$ точка u1 остается притягивающей, но значение производной в этой точке является предельным: | f'(u1) | = 1.

При значениях параметра логистического отображения $$\lambda = 1$$ и $$\lambda = 3$$ неподвижная точка этого отображения теряет устойчивость и появляется либо другая устойчивая неподвижная точка, как это произошло в первом случае, либо притягивающий цикл ; определение цикла будет дано ниже. Качественное изменение поведения решения (траектории отображения) при изменении параметра называется бифуркацией.

Пусть теперь $$3 < \lambda \le 1 + \sqrt{6} .$$ Как уже отмечалось, при значении параметра $$\lambda = 3$$ происходит бифуркация: неподвижная точка $$u_{2}= 1 - \lambda ^{ - 1}$$ из притягивающей превращается в отталкивающую: | f'(u) | > 1 при $$\lambda > 3.$$ После того как точка стала отталкивающей, рассмотрим корни u3, u4 уравнения f2(u) = u, или $$\lambda ^{2}u^{2} - \lambda (\lambda + 1)u + (\lambda + 1) = 0.$$

Заметим, что если u1 — предельная точка отображения f(u) = u, то она является также и предельной точкой отображения f2(u) = u. Действительно, f2(u1) = f(f(u1)) = f(u1) = u1, где u1 — любая предельная точка рассматриваемого отображения, отличная от корней уравнения f2(u) = u. Тогда, зная два корня уравнения f2(u) = u точки u3, u4 легко находятся как корни квадратного уравнения, они есть

Эти корни связаны соотношениями

f(u3) = u4, f(u4) = u3.

В данном случае говорят, что отображение имеет цикл периода 2, который будем обозначать P2. Его наличие, например, в популяционной модели говорит об изменении численности особей с периодом в 2 единицы времени. Траектория для случая такого цикла изображена на рис. 5.5. Можно считать, что неподвижная (предельная) точка отображения есть цикл периода 1.

Переход от цикла P1 (предельная точка логистического отображения) к циклу P2 называют бифуркацией рождения цикла (удвоения периода) .

Определение. Точка $$a \in X$$ называется периодической периода m, если fm (a) = a и $$f^{i} (a) \ne a$$ при 0 < i < m.

Отметим, что каковы бы ни были попарно различные точки u1, u2, ..., um, если положить f(ui) = ui + 1, i = 1, 2, ..., m - 1 и f(u_m) = u1, то рассматриваемое отображение будет иметь периодическую траекторию периода m: u1, u2, ..., um, u1, u2, ..., um, ...

Если к тому же f(u) имеет первую производную, то в окрестности каждой из точек ui выполнено

или

Будем рассматривать fm(u) как сложную функцию. Пользуясь правилом дифференцирования сложной функции, получим

Если $$\left|{\mathop \Pi\limits_{k = 1}^{m} f^{\prime}(u_k)}\right| < 1,$$ то траектория $$\{f_k (u_0 )\}_{k = 0}^\infty$$ приближается к циклу {u1, ..., uk}, или $$\{u^{k}\}_{k = 1}^{m} .$$ Такой цикл называется притягивающим циклом, а величина $$\left|{\mathop \Pi\limits_{k = 1}^{m} f^{\prime}(u_k)}\right|$$ — мультипликатором цикла. Цикл может быть как притягивающим, так и отталкивающим.

Определение. Цикл Pm = {u1, ..., um} отображения $$f: X \rightarrow X,$$ переводящего множество X в себя, называется притягивающим, если существует число k0, такое, что для любого k > k0 траектория $$\{f_k (u_0)\}_{k = 0}^\infty$$ распадается на m последовательностей, каждая из которых сходится к точкам u1, ..., um соответственно.

Достаточным условием существования притягивающего (отталкивающего) цикла является выполнение неравенства $$\mu (P_{m}) > 1,$$ где $$\mu (P_m) = {\mathop \Pi\limits_{k = 1}^{m} f^{\prime}(u_k)}, u \in P_m$$ — мультипликатор цикла.

Отметим интересные свойства функции f2(u), в частности, ее график пересекается с прямой y = u не только в неподвижных точках рассматриваемого отображения, u1, u2, но и в точках цикла P2. Таким образом, можно сказать, что бифуркация рождения цикла обусловлена потерей устойчивости одной предельной точки и появлением двух устойчивых предельных точек отображения f2(u). На рис. 5.6а, в показано поведение функции f2(u) при разных значениях параметра $$\lambda$$ ( $$\lambda = 2, 8; \lambda = 1 + \sqrt{5}$$ ).

При увеличении $$\lambda$$ у отображения появляются новые неподвижные точки. Мультипликатор цикла P2 вычисляется следующим образом:

Очевидно, что $$| \mu (u_{3}, u_{4}) | < 1,$$ если $$3 < \lambda < 1 + \sqrt{6},$$ тогда цикл P2 — притягивающий. Траектория $$\{f_k(u_0 )\}_{k = 0}^\infty$$ притягивается циклом {u3, u4} и подпоследовательность $$\{f^{2k} (u_0 )\}_{k = 0}^\infty$$ сходится к одной точке цикла, а $$\{f^{2k + 1} (u_0 )\}_{k = 0}^\infty$$ — к другой.

Знак мультипликатора дает информацию о характере приближения траектории к циклу. В частности, если $$3 < \lambda < 1 + \sqrt{5} ,$$ то подпоследовательности $$\{f^{2k} (u_0 )\}_{k = 0}^\infty $$ и $$\{f^{2k + 1} (u_0 )\}_{k = 0}^\infty,$$ начиная с некоторого u, являются монотонными, одна из них возрастающая, а другая — убывающая, что зависит от знаков f'(u3) и f'(u4).

При $$1 + \sqrt{5} < \lambda \le 1 + \sqrt{6}$$ значение мультипликатора $$\mu < 0,$$ и подпоследовательности $$\{f^{2k} (u_0 )\}_{k = 0}^\infty $$ и $$\{f^{2k + 1} (u_0 )\}_{k = 0}^\infty$$ приближаются к точкам {u3, u4} немонотонно.

Рассмотрим теперь случай $$1 + \sqrt{6} \le \lambda < 3, 54 \ldots$$

При $$\lambda = 1 + \sqrt{6}$$ происходит вторая бифуркация удвоения периода.

Цикл {u3, u4} из притягивающего превращается в отталкивающий, $$| \mu (u_{3}, u_{4}) | > 1$$ при $$\lambda > 1 + \sqrt{6} .$$ Появляется новый притягивающий цикл P4:

Для популяционной динамики это означает, что численность особей колеблется с периодом 4 единицы времени. Соответствующий график приведен на рис. 5.7.

При $$\lambda \approx 3, 54$$ цикл P4 периода 4 становится отталкивающим, $$| \mu (u_{5}, \dots , u_{8}) | > 1$$ ; при этом появляются притягивающий цикл $$P_8$$ периода 8. Дальнейшее увеличение параметра $$\lambda$$ будет приводить к появлению циклов P16, P32 и т.д. Происходит каскад бифуркаций удвоения периода.

Заметим, что рассмотренный простой процесс имеет сложное поведение. Наблюдается каскад бифуркаций при увеличении величины $$\lambda$$ ; кроме того, все циклы, которые при этом встречаются, имеют период 2p. Эта важнейшая закономерность, которая прослеживается не только в расчетах, но и в природе! Рассмотренные бифуркации при увеличении $$\lambda$$ можно наглядно представить на бифуркационной диаграмме ( рис. 5.8). Диаграмма получается, если обозначить через $$\Lambda _{1}, \Lambda _{2}$$ те значения $$\lambda,$$ в которых происходят бифуркации, а через $$\lambda _{1}, \lambda _{2}, ...$$ при которых u = 0, 5 является элементом циклов P2, P4, ...; по вертикальной оси откладываются значения предельных точек отображения. Обозначим за d1, d2, ... величины, равные расстоянию между x = 0, 5 и ближайшим к нему элементом цикла P2 при $$\lambda = \lambda _{k}.$$ Численный эксперимент показал, что $$\Lambda _{k}$$ и $$\lambda _{k}$$ при достаточно больших k ведут себя, как геометрическая прогрессия со знаменателем $$\delta = 4, 66920\ldots,$$ т.е.

Отношение dk/dk + 1 имеет предел, равный $$\alpha = 2, 50290 \ldots .$$ Эти закономерности были замечены американским математиком Фейгенбаумом.

При дальнейшем увеличении $$\lambda$$ последовательность $$\left\{{u_k}\right\}_{k = 0}^\infty$$ приобретает хаотический характер ( $$\lambda = \lambda_\infty \approx 3, 569$$ ), что видно на рис. 5.9.

Примечательно, что каскады Фейгенбаума имеют фрактальный характер (т.е. сохраняют подобие при изменении масштабов, рис. 5.10 а, б).

Изучение графиков функций f2(u) и f1(u) показывает, что их фрагменты вблизи максимумов близки друг к другу, более того, они отличаются лишь масштабами. Оказывается, что такое же подобие имеет место для функции $$f^{2^k}, n > 1$$ при $$\lambda = \lambda _{k},$$ и выполняется тем точнее, чем больше n. Если положить u' = u - 1 / 2 (в дальнейшем штрих будем опускать) и считать $$\alpha$$ коэффициентом растяжения вдоль осей, то для некой симметричной функции g(u), определенной на отрезке [- 1, 1], можно получить следующее функциональное уравнение:

которое универсально определяет $$\alpha:$$

Вблизи максимума g(x) должна быть близка к квадратичной параболе, причем g(0) = 1. В теории универсальности показывается, что эта функция вычисляется с помощью ряда

g(u) = 1 - 1, 52763u2 + 0, 104815u4 - 0, 0267057u6 + ...

Пусть теперь $$\lambda = 3, 83.$$ В этом случае из хаотической области, изображенной на рис. 5.10, появляется устойчивый цикл P3 (рис. 5.11 а, b представляют циклы в последовательные моменты времени).

Циклу на рисунке выше соответствует самое большое окно устойчивых циклов $$P_{3 \cdot 2^{k}}.$$ Чередование хаотических и регулярных зон — называется перемежаемостью. Возможно, нечто подобное наблюдается в гидродинамических потоках, где ламинарные зоны чередуются с турбулентными.

5.6. Задачи

Методами простых итераций и Ньютона решить уравнение

$$x^2 - e^{- x} = 0, x_0 \in \left[{0.5, 1}\right].$$

Решение. Метод простых итераций будет $$$ x_{k + 1} = e^{- \frac{x_k}{2}}, x_0 = 0, 75 $,$$ $$$ F(x) = e^{- \frac{x_k}{2}}, \left|{F^{\prime}_x (x)}\right| < 1 $$$ при $$x \in \left[{0.5, 1}\right].$$

Приведем таблицу приближений до точности | x5 - x4 | < 10 - 4.

k	0	1	2	3	4	5
xk	0, 75	0, 6873	0, 7091	0, 7015	0, 7042	0, 7032

Метод Ньютона запишется как

$$$ x_{k + 1} = x_k - \frac{{(x_k)}^2 - e^{- x_k}}{2x_k - e^{- x_k}}, x_0 = 1. $$$

Результаты расчетов:

k	0	1	2	3
xk	1	0, 7330	0, 7038	0, 7035

Построить итерационные методы для вычисления корней уравнения

x3 + 3x2 - 1 = 0.

Использовать метод простых итераций.

Решение. Все три корня данного уравнения лежат на отрезках [- 3, - 2], [- 1, 0], [0, 1].

Построим итерационный процесс для вычисления первого корня, лежащего на отрезке [- 3, - 2]:

$$$ x_{k + 1} = x_k ^{- 2} - 3, F(x) = x^{- 2} - 3, \left|{F^{\prime}_x (x)}\right| = \left|{- 2x^{- 3}}\right| \le \frac{1}{4} < 1, $$$

т.е. для начального приближения $$x_0 \in \left[{- 3, - 2}\right]$$ метод простых итераций сходится.

Для вычисления двух оставшихся корней, лежащих на отрезках [- 1, 0] и [0, 1], построим итерационный метод

$$x_{k + 1} = \pm (x_k + 3)^{- \frac{1}{2}}, F(x) = \pm (x_k + 3)^{- \frac{1}{2}}. $$$

Поскольку для рассматриваемых отрезков $${\left|{F^{\prime}_x (x)}\right| = \left|{2(x + 3)^{- \frac{1}{2}}}\right|^3 < 1} $,$$ то этот итерационный процесс сходится к корням рассматриваемого уравнения.

Исследовать сходимость итерационного метода релаксации для численного решения уравнения вида f(x) = 0:$$x_{k + 1} = x_{k} - \tau f(x_{k}), x_{0} = a, f'(x) > 0.$$

Решение. Итерационный метод релаксации записывается как

xk + 1 = F(xk), x0 = 0,

$$F(x) = x - \tau f(x).$$

Из условия того, что отображение является сжимающим, получим $$\left|{F(x) - F(y)}\right| = \left|{x - y - \tau (f(x) - f(y))}\right| \le q\left|{x - y}\right|,$$ где $$q = \max \left|{1 - \tau f^{\prime}(x + \alpha (x - y))}\right|,$$ $$0 \le \alpha \le 1.$$

Положим $$0 < f^{\prime}_{\min } \le f^{\prime} \le f^{\prime}_{\max },$$ тогда $${q}(\tau ) = \max \left\{\left|{1 - \tau f^{\prime}_{\min }}\right|\right.,$$ $$\left.\left|{1 - \tau f^{\prime}_{\max}}\right|\right\}.$$

Отсюда видно, что итерационный процесс будет сходиться, если $$\left|{1 - \tau f^{\prime}_{\max }}\right| < 1,$$ или $$\tau < 2/f^{\prime}_{\max } ,$$ а оптимальное значение параметра, при котором q = qmin достигается при $$\left|{1 - \tau f^{\prime}_{\min }}\right| = \left|{1 - \tau f^{\prime}_{\max }}\right|,$$ или $$\tau = \tau_0 = 2\left({f^{\prime}_{\min } + f^{\prime}_{\max }}\right)^{- 1}, q(\tau_0) = (1 - f^{\prime}_{\min } /f^{\prime}_{\max })(1 + f^{\prime}_{\min } /f^{\prime}_{\max }).$$

Построить метод простых итераций и итерационный процесс Ньютона для системы уравнений

f(x, y) = x + 3 lg x - y2 = 0, 
 g(x, y) = 2x2 - xy - 5x + 1 = 0.

Решение. Графики функций f(x, y) и g(x, y) приведены на рис. 5.12.

(рис 5.12)

Метод простых итераций запишем как

$$\begin{gather*} x_{k + 1} = y_k ^2 - 3\lg x_k, \\ y_{k + 1} = 2x_k + \frac{1}{x_k} - 5, \end{gather*}$$

где $$$ x_0 = 3, 4, y_0 = 2, 2, F(x, y) = y^2 - 3 \lg x; G(x, y) = 2x + \frac{1}{x} - 5 $.$$

Матрица Якоби для такого процесса запишется как

$$$ {\mathbf{J}} = \left( \begin{array}{cc} {\frac{\partial F}{\partial x}} {\frac{\partial F}{\partial y}} \\ {\frac{\partial G}{\partial x}} {\frac{\partial G}{\partial y}} \\ \end{array} \right) = \left( \begin{array}{cc} {- \frac{3\lg e}{x}} {2y} \\ {2 - x^{- 2}} 0 \\ \end{array} \right). $$$

Несложно проверить, вычислив норму матрицы Якоби, что для приведенного начального приближения достаточное условие сходимости не выполняется.

Рассмотрим другой итерационный процесс:

$$\begin{gather*} x_{k + 1} = \left({\frac{x_k (y_k + 5) - 1}{2}}\right)^{1/2}, \\ y_{k + 1} = \left({x_k + 3\lg x_k}\right)^{1/2}, \\ F(x, y) = \left({\frac{x_k (y_k + 5) - 1}{2}}\right)^{1/2}, G(x, y) = \left({x + 3\lg x}\right)^{1/2}. \end{gather*}$$

Матрица Якоби для этого итерационного метода будет

$$$ {\mathbf{J}} = \left( \begin{array}{cc} {\frac{5 + y}{2\sqrt{2} \sqrt{x(y + 5) - 1}}} {\frac{x}{{2\sqrt{2} \sqrt {x(y + 5) - 1}}}} \\ {\frac{{1 + \frac{{3\lg e}}{x}}}{{2\sqrt {x + 3\lg x}}}} 0 \\ \end{array} \right). $$$

В окрестности начального приближения условие сходимости выполнено. Таблица первых пяти приближений будет

k	0	1	2	3	4	5
x	3, 4	3, 426	3, 451	3, 466	3, 475	3, 480
y	2, 2	2, 243	2, 2505	2, 255	2, 258	2, 259

Построить итерационный метод Ньютона для вычисления $$\sqrt[n]{a},$$ a > 0, $$n \in R.$$

Решение. Найдем корень уравнения

f(x) = xn - a = 0.

Метод Ньютона для этого уравнения запишется

$$\begin{gather*} x_{k + 1} = x_k - \frac{{f\left({x_k}\right)}}{{f^{\prime}\left({x_k}\right)}} = x_k - \frac{{x_n^{k} - a}}{{nx_m^{p - 1}}} = \\ = \frac{{n - 1}}{n}x_k + \frac{a}{{nx_k^{n - 1}}} = \frac{1}{n}\left[{\left({n - 1}\right)x_k + \frac{a}{{x_k^{n - 1}}}}\right]. \end{gather*}$$

В частности, при n = 2 имеем

$$$ x_{k + 1} = \frac{1}{2}\left({x_k + \frac{a}{x_k}}\right). $$$

Получить расчетные формулы метода Ньютона для численного решения системы двух нелинейных уравнений f(x, y) = 0, g(x, y) = 0.

Решение. Положим $$x_{k + 1} = \Delta x_{k},$$ $$y_{k + 1} = y_{k} + \Delta y\_ k,$$ f(xk, yk) = fk, получим

$$\begin{gather*} f(x_k , y_k) + \left({\frac{\partial f}{\partial x}}\right)^{k} \Delta x_k + \left({\frac{\partial f}{\partial y}}\right)^{k} \Delta y_k = 0, \\ g(x_k , y_k) + \left({\frac{\partial g}{\partial x}}\right)\Delta x_k + \left({\frac{\partial g}{\partial y}}\right)\Delta y_k = 0, \end{gather*}$$

откуда

$$\begin{gather*} \Delta x_k = \frac{{\left|{{\mathbf{X}}_k}\right|}}{{\left|{{\mathbf{J}}_k}\right|}} = \left| \begin{array}{cc} {\left({\frac{{\partial f}} {{\partial x}}}\right)_k} {\left({\frac{{\partial f}} {{\partial y}}}\right)_k} \\ {\left({\frac{{\partial g}} {{\partial x}}}\right)_k} {\left({\frac{{\partial g}} {{\partial y}}}\right)_k} \\ \end{array} \right|^{- 1}\cdot \\ \cdot \left| \begin{array}{cc} {- f_k} {\left({\frac{\partial f}{\partial y}}\right)_k} \\ {- g_k} {\left({\frac{\partial g}{\partial y}}\right)_k} \\ \end{array} \right| = \frac{{- f_k \left({\frac{\partial g}{\partial y}}\right)_k + g_k \left({\frac{\partial f}{\partial y}}\right)_k}} {{\left({\frac{\partial f}{\partial x}}\right)_k \left({\frac{\partial g}{\partial y}}\right)_k - \left({\frac{\partial g}{\partial x}}\right)_k \left({\frac{\partial f}{\partial y}}\right)_k}}, \\ \Delta y_k = \frac{{\left|{{\mathbf{Y}}_k}\right|}} {{\left|{{\mathbf{J}}_k}\right|}} = \left| \begin{array}{cc} {\left({\frac{{\partial f}} {{\partial x}}}\right)_k} {\left({\frac{{\partial f}} {{\partial y}}}\right)_k} \\ {\left({\frac{{\partial g}} {{\partial x}}}\right)_k} {\left({\frac{{\partial g}} {{\partial y}}}\right)_k} \\ \end{array} \right|^{- 1} \cdot \\ \cdot {\left| \begin{array}{cc} {\left({\frac{{\partial f}} {{\partial x}}}\right)_k} {- f_k} \\ {\left({\frac{{\partial g}} {{\partial x}}}\right)_k} {- g_k} \\ \end{array} \right|} = \frac{{- g_k \left({\frac{{\partial g}} {{\partial x}}}\right)_k + f_k \left({\frac{{\partial f}} {{\partial x}}}\right)_k}} {{\left({\frac{{\partial f}} {{\partial x}}}\right)_k \left({\frac{{\partial g}} {{\partial y}}}\right)_k - \left({\frac{{\partial g}} {{\partial x}}}\right)_k \left({\frac{{\partial f}} {{\partial y}}}\right)_k}}, \end{gather*}$$

где $$$ \mathbf{J}_k = \left( \begin{array}{cc} {\left({\frac{{\partial f}} {{\partial x}}}\right)_k} {\left({\frac{{\partial f}} {{\partial y}}}\right)_k} \\ {\left({\frac{{\partial g}} {{\partial x}}}\right)_k} {\left({\frac{{\partial g}} {{\partial y}}}\right)_k} \\ \end{array} \right) $$$ — матрица Якоби, $${\mathbf{X}}_k , {\mathbf{Y}}_k$$ — матрицы

$$$ {\mathbf{X}}_k = \left( \begin{array}{cc} {- f_k} {\left({\frac{{\partial f}} {{\partial x}}}\right)_k} \\ {- g_k} {\left({\frac{{\partial f}} {{\partial y}}}\right)_k} \\ \end{array} \right), {\mathbf{Y}}_k = \left( \begin{array}{cc} {\left({\frac{{\partial f}} {{\partial x}}}\right)_k} {- f_k} \\ {\left({\frac{{\partial f}} {{\partial y}}}\right)_k} {- g_k} \\ \end{array} \right). $$$

Тогда запишем расчетные формулы для итерационного метода Ньютона:

$$x_{k + 1} = x_{k} + \Delta x_{k} y_{k + 1} = y_{k} + \Delta y_{k}.$$

Найти приближенное решение системы двух нелинейных алгебраических уравнений x3 - y3 - 1 = 0, xy3 - y - 4 = 0, используя метод Ньютона.

Решение.

Пусть f(x, y) = x3 - y3 - 1, g(x, y) = xy3 - y - 4. Начальное приближение можно найти графически: x0 = y0 = 1, 5.

Матрицы $$\mathbf{J}_k, \mathbf{X}_k, \mathbf{Y}_k$$ в этом случае будут

$${\mathbf{J}}_k = - \left( \begin{array}{cc} {3x_k^2} {- 2y_k} \\ {y_k^3} {3x_k y_k^2 - 1} \\ \end{array} \right), {\mathbf{X}}_k = - \left( \begin{array}{cc} {x_k^3 - y_k^3 - 1} {- 2y_k} \\ {x_k y_k^3 - y_k - 4} {3x_k y_k^2 - 1} \\ \end{array} \right), \\ {\mathbf{Y}}_k = - \left( \begin{array}{cc} {3x_k^2} {x_k^3 - y_k^3 - 1} \\ {y_n^3} {x_n y_n^3 - y_n - 4} \\ \end{array} \right).$$

Следующее приближение вычисляется по формуле Ньютона:

$$$ x_{k + 1} = x_k + \frac{{\left|{{\mathbf{X}}_k}\right|}} {{\left|{{\mathbf{J}}_k}\right|}}, y_{k + 1} = y_k + \frac{{\left|{{\mathbf{Y}}_k}\right|}} {{\left|{{\mathbf{J}}_k}\right|}}. $$$

Результаты вычислений первых двух итераций приведены в таблице (точность $$\varepsilon \approx 10^{ - 3}$$ ) .

k	xk; yk	fk; gk	$${\left|\mathbf{J}_k\right|}$$	$${\left|\mathbf{X}_k\right|}$$	$${\left|\mathbf{Y}_k\right|}$$
0	1,5; 1,5;	0,12500 -4,33750	71,71875	-0,171875	-3,3750
1	1,502397 1,547059	-0,002170 0,015844	77,73277	0,0277988	0,1153255
2	1,5020396 1,545570	0,0000017 0,000019

Найти приближенное решение системы трех нелинейных уравнений x2 + y2 + z2 = 1, 2x2 + y2 - 4z = 0, 3x2 - 4y + z2 = 0, используя метод Ньютона. За начальное приближение решения выбрать точку x0 = y0 = z0 = 0, 5.

Решение. Обозначим $${\mathbf{F}} = (x^2 + y^2 + z^2 - 1, 2x^2 + y^2 - 4z, 3x^2 - 4y + z^2 )^T ,$$ тогда матрица Якоби будет $${\mathbf{J}} = \left( \begin{array}{ccc} {2x} {2y} {2z} \\ {4x} {2y} {- 4} \\ {6x} {- 4} {2z} \end{array} \right).$$ Соответственно, для точки начального приближения вычислим

$$$ {\mathbf{F}}_0 = \left( \begin{array}{l} 0, 25 + 0, 25 + 0, 25 - 1 \\ 0, 5 + 0, 25 - 2, 00 \\ 0, 75 - 2, 00 + 0, 25 \end{array} \right) = \left( \begin{array}{l} - 0, 25 \\ - 1, 25 \\ - 1, 00 \\ \end{array} \right), \\ {\mathbf{J}}_0 = \left( \begin{array}{ccc} 1 1 1 \\ 2 1 {- 4} \\ 3 {- 4} 1 \\ \end{array} \right), \det {\mathbf{J}}_0 = - 40 \ne 0, \\ {\mathbf{J}}_0^{ - 1} = - \frac{1} {{40}}\left( \begin{array}{ccc} {- 15} {- 5} {- 5} \\ {- 14} {- 2} 6 \\ {- 11} 7 {- 1} \\ \end{array} \right). $$$

Вычислим первое приближение

$$$ {\mathbf{x}}_1 = {\mathbf{x}}_0 - {\mathbf{J}}_0^{- 1}{\mathbf{F}}_0 = \left( \begin{array}{l} 0, 5 \\ 0, 5 \\ 0, 5 \\ \end{array} \right) - \\ - \frac{1} {{40}}\left( \begin{array}{ccc} {- 15} {- 5} {- 5} \\ {- 14} {- 2} 6 \\ {- 11} 7 {- 1} \\ \end{array} \right) \left( \begin{array}{l} - 0, 25 \\ - 1, 25 \\ - 1, 00 \\ \end{array} \right) = \left( \begin{array}{l} 0, 875 \\ 0, 500 \\ 0375 \\ \end{array} \right). $$$

Для второго приближения

$$$ {\mathbf{F}}_1 = \left( \begin{array}{l} 0, 875^2 + 0, 500^2 + 0, 375^2 - 1 \\ 2 \cdot 0, 875^2 + 0, 500^2 - 4 \cdot 0, 375 \\ 3 \cdot 0, 875^2 - 4 \cdot 0, 500 + 0, 375^2 \\ \end{array} \right) = \left( \begin{array}{l} 0, 15625 \\ 0, 28125 \\ 0, 43750 \\ \end{array} \right), \\ {\mathbf{J}}_1 = \left( \begin{array}{ccc} {2 \cdot 0, 875} {2 \cdot 0, 500} {2 \cdot 0, 375} \\ {4 \cdot 0, 875} {2 \cdot 0, 500} {- 4} \\ {6 \cdot 0, 875} {- 4} {2 \cdot 0, 375} \\ \end{array} \right) = \left( \begin{array}{ccc} {1, 750} 1 {0, 750} \\ {3, 500} 1 {- 4} \\ {5, 250} {- 4} {0, 75} \\ \end{array} \right), \\ \det {\mathbf{J}}_1 = - 64, 75, \\{\mathbf{J}}_1^{- 1} = - \frac{1} {{64, 75}}\left( \begin{array}{ccc} {- 15, 25} {- 3, 75} {- 4, 75} \\ {- 23, 625} {- 2, 625} {9, 625} \\ {- 19, 25} {12, 25} {- 1, 75} \\ \end{array} \right). $$$

По формуле Ньютона получим

$$$ {\mathbf{x}}_2 = {\mathbf{x}}_1 - {\mathbf{J}}_1^{- 1} {\mathbf{F}}_1 = \left( \begin{array}{l} 0, 875 \\ 0, 500 \\ 0, 375 \\ \end{array} \right) + \\ + \frac{1} {{64, 75}}\left( \begin{array}{ccc} {- 15, 25} {- 3, 75} {- 4, 75} \\ {- 23, 625} {- 2, 625} {9, 625} \\ {- 19, 25} {12, 25} {- 1, 75} \\ \end{array} \right) \cdot \left( \begin{array}{l} 0, 15625 \\ 0, 28125 \\ 0, 43750 \\ \end{array} \right) = \left( \begin{array}{l} 0, 78981 \\ 0, 49662 \\ 0, 36993 \\ \end{array} \right). $$$

Аналогично находим третье приближение, которым и ограничимся:

$${\mathbf{x}}_3 = \left( \begin{array}{l} 0, 78521 \\ 0, 49662 \\ 0, 36992 \\ \end{array} \right).$$

5.7. Задачи для самостоятельного решения

x - 2x - 1 = 0,