Рассмотрим систему из двух уравнений $F_1(x,y) =0,\;F_2(x,y)=0$ , обладающую свойством покоординатной монотонности: с ростом функции также растут. Задача состоит в нахождении всех корней системы.

Введение

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

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

Пример

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

$F_1(x,y)=x+y+\varphi(x)-0.1,\\F_2(x,y)=x+y-\varphi(x)+0.1,\\\varphi(x)=e^{-x^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).

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

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

Алгоритм

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

$\widehat M_{ji} \leq \frac{\partial F_j(\mathbf x)}{\partial x_i}\leq \widetilde M_{ji},$

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

Пусть $\mathbf x^{(k)}$ - текущее приближение, в котором функции и имеют разные знаки (на рисунке такая точка будет находиться между двумя кривыми). Тогда, чтобы, во-первых, зафиксировать направление движения алгоритма (влево-вверх или вправо-вниз) и, во-вторых, гарантировать, что следующее приближение $\mathbf x^{(k+1)}$ останется в той же области и не перескочет корень, оценим приращение значений функций по формуле конечных приращений. Получим

$F_1(\mathbf x^{(k+1)})=F_1(\mathbf x^{(k)})+\frac{\partial F_1(\widetilde {\mathbf x})}{\partial x}(x^{(k+1)}-x^{(k)})+\frac{\partial F_1(\widetilde {\mathbf x})}{\partial y}(y^{(k+1)}-y^{(k)}),\\F_2(\mathbf x^{(k+1)})=F_2(\mathbf x^{(k)})+\frac{\partial F_2(\widehat {\mathbf x})}{\partial x}(x^{(k+1)}-x^{(k)})+\frac{\partial F_2(\widehat {\mathbf x})}{\partial y}(y^{(k+1)}-y^{(k)}).$

Тогда, если мы хотим, чтобы алгоритм двигался влево-вверх, то есть $x^{(k+1)}<x^{(k)}, \; y^{(k+1)}>y^{(k)}$ , то для сохранения свойства нужно оценить элементы якобиана снизу и сверху:

$F_1(\mathbf x^{(k+1)})\geq F_1(\mathbf x^{(k)})+\widetilde M_{11}(x^{(k+1)}-x^{(k)})+\widehat M_{12} (y^{(k+1)}-y^{(k)}),\\F_2(\mathbf x^{(k+1)})\leq F_2(\mathbf x^{(k)})+\widehat M_{21}(x^{(k+1)}-x^{(k)})+\widetilde M_{22}(y^{(k+1)}-y^{(k)}).$

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

$\widetilde M_{11}(x^{(k+1)}-x^{(k)})+\widehat M_{12} (y^{(k+1)}-y^{(k)}) = -F_1(\mathbf x^{(k)}),\\\widehat M_{21}(x^{(k+1)}-x^{(k)})+\widetilde M_{22}(y^{(k+1)}-y^{(k)}) = -F_2(\mathbf x^{(k)})$

относительно $x^{(k+1)},y^{(k+1)}$ , то получим, что

$F_1(\mathbf x^{(k+1)})\geq 0, \; F_2(\mathbf x^{(k+1})\leq 0.$

Но для этого нужно, чтобы решение данной линейной системы соответствовало выбранному направлению $x^{(k+1)}<x^{(k)},\; y^{(k+1)}>y^{(k)}$ . Это будет достигнуто при определенной ориентации векторов, а именно когда

$\left|\matrix{\widetilde M_{11} & \widehat M_{12}\\ \widehat M_{21} & \widetilde M_{22}}\right| > 0.$

Аналогично в случае $x^{(k+1)}>x^{(k)},\; y^{(k+1)}<y^{(k)}$ , решая систему:

$\widehat M_{11}(x^{(k+1)}-x^{(k)})+\widetilde M_{12} (y^{(k+1)}-y^{(k)}) = -F_1(\mathbf x^{(k)}),\\\widetilde M_{21}(x^{(k+1)}-x^{(k)})+\widehat M_{22}(y^{(k+1)}-y^{(k)}) = -F_2(\mathbf x^{(k)})$

относительно $x^{(k+1)},y^{(k+1)}$ , получаем, что

$F_1(\mathbf x^{(k+1)})\leq 0, \; F_2(\mathbf x^{(k+1})\geq 0.$

поскольку

$F_1(\mathbf x^{(k+1)})\leq F_1(\mathbf x^{(k)})+\widehat M_{11}(x^{(k+1)}-x^{(k)})+\widetilde M_{12} (y^{(k+1)}-y^{(k)}),\\F_2(\mathbf x^{(k+1)})\geq F_2(\mathbf x^{(k)})+\widetilde M_{21}(x^{(k+1)}-x^{(k)})+\widehat M_{22}(y^{(k+1)}-y^{(k)}).$

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

$\left|\matrix{\widehat M_{11} & \widetilde M_{12}\\ \widetilde M_{21} & \widehat M_{22}}\right| < 0.$

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

Код для вычислений

# оценки якобиана
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 итераций.

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

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