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

Введение

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

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

Пример

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

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 итераций.

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