Как найти все корни системы 2×2?

Рассмотрим систему из двух уравнений , обладающую свойством покоординатной монотонности: с ростом функции также растут. Задача состоит в нахождении всех корней системы.

Введение

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

В настоящей статье предложена разновидность метода Ньютона, которая позволяет зафиксировать направление движения алгоритма с помощью использования подходящих оценок якобиана системы уравнений. С помощью этой техники сканируется всё интересующее пространство поиска. По сравнению с методом «лесенки», в котором двумерная задача сводится к одномерной бисекции по чередующимся координатам, метод позволяет ускорить сходимость.

Пример

Рассмотрим следующую систему:

function ret = phi(x)   ret = exp(-x^2); endfunction  function ret = dphi(x)   ret = -2 * x * exp(-x^2); endfunction  function ret = F1(x, y)   ret = x + y + phi(x) - 0.1; endfunction  function ret = F2(x, y)   ret = x + y - phi(x) + 0.1; endfunction  X = linspace(-3, 3, 100); Y = linspace(-3, 3, 100); [X, Y] = meshgrid(X, Y); Z1 = arrayfun(@F1, X, Y); Z2 = arrayfun(@F2, X, Y); hold on contour(X, Y, Z1, [0 0]) contour(X, Y, Z2, [0 0])

function [z, jac] = f(p)   x = p(1);   y = p(2);    z(1) = F1(x, y);   z(2) = F2(x, y);   jac(1, 1) = 1 + dphi(x);   jac(1, 2) = 1;   jac(2, 1) = 1 - dphi(x);   jac(2, 2) = 1; endfunction  x_init = [-0.5; 0.5]; [sol, fval, info] = fsolve (@f, x_init, optimset ("jacobian", "on")) x_init = [0.5; 0.5]; [sol, fval, info] = fsolve (@f, x_init, optimset ("jacobian", "on")) x_init = [0; 0]; [sol, fval, info] = fsolve (@f, x_init, optimset ("jacobian", "on")) 

Запустим стандартный решатель fsolve (среда Octave).

При начальном приближении получаем корень .

При начальном приближении получаем корень .

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

Алгоритм

Допустим, что нам известны оценки снизу и сверху для элементов якобиана системы:

причем неотрицательные матрицы обладают нужными нам свойствами (об этом дальше).

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

Тогда, если мы хотим, чтобы алгоритм двигался влево-вверх, то есть

Значит, если мы решим такую систему:

относительно , то получим, что

Но для этого нужно, чтобы решение данной линейной системы соответствовало выбранному направлению

относительно , получаем, что

поскольку

При этом должно быть выполнено условие

Запустим алгоритм для нашего примера.

# оценки якобиана J_upper = [2, 1; 2, 1]; J_lower = [0, 1; 0, 1]; Q1 = [J_upper(1, 1), J_lower(1, 2); J_lower(2, 1), J_upper(2, 2)]; Q2 = [J_lower(1, 1), J_upper(1, 2); J_upper(2, 1), J_lower(2, 2)];  # алгоритм, который идет влево-вверх disp("Algorithm 1") num_iter = 30; x = x_init; for i = 1:num_iter   x -= Q1 ^ -1 * [F1(x(1), x(2)); F2(x(1), x(2))] endfor  # алгоритм, который идет вправо-вниз disp("Algorithm 2") x = x_init; for i = 1:num_iter   x -= Q2 ^ -1 * [F1(x(1), x(2)); F2(x(1), x(2))] endfor

Получим в одном случае корень , в другом случае корень за 30 итераций.

Таким образом, предложенный алгоритм позволяет избавиться от многократного применения сравнительно трудоемкого одномерного метода бисекции в алгоритме «лесенка».