Часть 2. Алгоритм Тарьяна для приведения нелинейной системы уравнений к вычисляемой последовательности подстановок

Это вторая часть статьи, посвященную вопросу применения алгоритма Тарьяна для решения систем уравнений. В первой части рассматривалась проблема поиска минимального набора уравнений, из неё нам понадобится только несколько определений.

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

Формализация задачи и вводный пример

Введём новое определение: вычисляемая последовательность формальных подстановок — каждый элемент последовательности (уравнение) содержит в правой части только уже вычисленные неизвестные, т.е. выраженные через известные переменные.

Зададим матрицу $S_{n \times n}$ для правой части по принципу из предыдущей статьи: $s_{ij} = 1$ , если -переменная входит в -ое уравнение.

Если последовательность вычисляема, то матрица — нижнетреугольная с нулями на диагонали. Поясним на примере:

$\begin{cases} x_1= A_1 \\ x_2 = \cos{x_1} + A_2 \\ x_3 = x_1^2 + x_1 + x_2 \\ x_4 = A_3 \end{cases}$

Соответствующая матрица :

$S = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{pmatrix}$

Декларируемая вычислимость состоит в том, что известный x_1 можно подставить в следующее за ним уравнение, найдя x_2 , известные x_1, x_2 подставляются в уравнение для x_3 и т.д.

Подстановки можно выразить в виде уравнения x = Sx + C . Решая уравнение (например, методом Гаусса), мы установим являются ли подстановки вычислимой последовательностью.

Алгоритм

Свести каждое уравнение системы к виду $x_i = f(x_1, \dots, x_n, t)$ . не должен входить в правую часть выражения. Иначе, отнести это выражение в группу ‘уравнений’.
Найти компоненты сильной связности с помощью алгоритма Тарьяна.
Удалить алгебраические петли. Запустить заново алгоритм Тарьяна.
Компоненты размера $1 \times 1$ являются подстановками, остальные — уравнениями.
Упорядочить и получить формальную вычислимую последовательность.

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

Пример

Этот пример и идея алгоритма были взяты из статьи About detection substitutions in nonlinear algebraic equations with help of Tarjan’s algorithm. Isakov A.A., Senichenkov Yu.B.

Система:

$\left\{\begin{array}{l} x_{1}=-7x_{7}+x_{8}-2 \\ x_{2}=3x_{4} - 2x_{6}+1 \\ x_{3}=\frac{x_{5}}{2}+8 \\ x_{4}=\frac{x_{6}}{x_{1}}+12 \\ x_{5}=(x_{7}+x_{8}+x_{3})^2-23 \\ x_{6}=x_{3}+x_{2}-3 \\ x_{7}=2 x_{3}+2 \\ x_{8}=(x_{7}-2)^3+2 \end{array}\right.$

Визуализация графа связей между уравнениями и переменными:

Об удалении алгебраических петель

Алгебраическая петля это простой цикл в нашем графе, который находится методом библиотеки NetworkX nx.simple_cycles. Фактически, это когда в результате подстановки выясняется, что переменная зависит от себя, нарушая принцип, озвученный в п. 1. алгоритма.

$1. \ x_3 \Longrightarrow x_5 \Longrightarrow x_7 \Longrightarrow x_8 \Longrightarrow x_3$ $2. \ x_2 \Longrightarrow x_4 \Longrightarrow x_6 \Longrightarrow x_2$

Определение циклового порядка (, cycling order) := {количество входящих ребер из вершин компоненты сильной связности (КСС); количество исходящих ребер в вершины текущей КСС}.

Необходимо выбрать вершину с максимальный цикловым порядком, отнести её в разряд ‘уравнений’ и удалить из графа. CO(x_5) = 3 , равный максимальному, её и удалим.
Если вершин с максимальным цикловым порядком несколько, то выберем ту, что первая встретилась в алгоритме Тарьяна. В данном случае, удалим вершину x_2 с порядком 2, которая встретилась первой.

Результаты работы алгоритма

Передём к программированию алгоритма. Ввод данных:

M = np.array([     [0, 0, 0, 0, 0, 0, 1, 1],     [0, 0, 0, 1, 0, 1, 0, 0],     [0, 0, 0, 0, 1, 0, 0, 0],     [1, 0, 0, 0, 0, 1, 0, 0],     [0, 0, 1, 0, 0, 0, 1, 1],     [0, 1, 1, 0, 0, 0, 0, 0],     [0, 0, 1, 0, 0, 0, 0, 0],     [0, 0, 0, 0, 0, 0, 1, 0], ])  tsk = SubstitutionDetector(dim=8) tsk.set_matrix(M) tsk.break_loops() tsk.get_answer()

Вывод программы:

Для цикла [2, 4, 7, 6] будет удалена вершина: 4. Причина: Максимальный цикловый порядок Для цикла [1, 3, 5] будет удалена вершина: 1. Причина: Максимальный цикловый порядок, первым встретился в алгоритме Тарьяна Индексы уравнений: [4, 1] Индексы(!) подстановок: [2, 6, 7, 0, 5, 3] Индексы(!) подстановок, которые действительно есть в системе: [2, 6, 7, 0, 5, 3]

Следовательно, «уравнениями» будут выражения, в левой части которых стоит x_5, x_2 , а «подстановками» будут и помним, что порядок имеет значение.

Уравнения:

$\left\{\begin{array}{l}x_{2}=3x_{4} - 2x_{6}+1 \\x_{5}=(x_{7}+x_{8}+x_{3})^2-23 \\\end{array}\right.$

Подстановки:

$\left\{\begin{array}{l}x_{3}=\frac{x_{5}}{2}+8 \\x_{7}=2 x_{3}+2 \\x_{8}=(x_{7}-2)^3+2 \\x_{1}=-7x_{7}+x_{8}-2 \\x_{6}=x_{3}+x_{2}-3 \\x_{4}=\frac{x_{6}}{x_{1}}+12 \\\end{array}\right.$

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

$\left\{\begin{array}{l} x_{3}=f(x_5) \\ x_{7}=f_1(x_5) \\ x_{8}=f_2(x_5) \\ x_{1}=f_3(x_5) \\ x_{6}=g(x_2, x_5) \\ x_{4}=g_1(x_2, x_5) \\ \end{array}\right.$

Теперь можно решить первую подсистему, которая была отнесена к категории «уравнений». Когда x_2, x_5 будут найдены, то через них можно найти и оставшиеся неизвестные.

Код

Утилитные функции:

Hidden text

def tarjan(graph: dict[int, set[int]]):     """     Пример представления графа.     graph = {     0: {1},     1: {3},     2: {1, 4},     3: {0},     4: {2, 5},     5: set()     }     :param graph:     :return:     """     n = len(graph)     sccs = []     index = [0]     indexes = [-1] * n     lows = [float('Inf')] * n     S = []      def strongconnect(v):         indexes[v] = index[0]         lows[v] = index[0]         index[0] += 1         S.append(v)         for chld in graph[v]:             if indexes[chld] == -1:                 strongconnect(chld)                 lows[v] = min(lows[v], lows[chld])             elif chld in S:                 lows[v] = min(lows[v], lows[chld])          if lows[v] == indexes[v]:             scc = [v]             w = S.pop()             while w != v:                 if w not in scc:                     scc.append(w)                 w = S.pop()             sccs.append(scc)      for v in graph.keys():         if indexes[v] == -1:             strongconnect(v)      return sccs    def matr2dct(m):     gr = {}     num_nodes = m.shape[0]     for i in range(num_nodes):         neighbors = set(np.nonzero(m[i])[0])         gr[i] = neighbors     return gr

Класс SubstitutionDetector:

Hidden text

import networkx as nx import numpy as np from tarjan.utils import tarjan, matr2dct   class SubstitutionDetector:     def __init__(self, dim):         self.translator = dict()         self.inverted_translator = dict()         self.graph = None         self.scc = None         self.matrix = np.zeros((dim, dim))         self.equations_indices = []         self.substitution_indices = []         self.in_system = []      def add_equation(self, varnum: int, dependnums: list[int]):         for i in dependnums:             self.matrix[varnum - 1][i - 1] = 1             self.in_system.append(varnum - 1)      def set_matrix(self, matr):         self.matrix = matr         for i in range(self.matrix.shape[0]):             if any(elem != 0 for elem in self.matrix[i]):                 self.in_system.append(i)         self.initialize()      def cycling_order(self, vertex_idx):         """         max(исходящие ребра, выходящие ребра) в вершину той же компоненты сильной связности         :param vertex_idx:         :return:         """         group = self.get_component_by_vertex(vertex_idx)         ones_in_row = 0         ones_in_column = 0          for i, r in enumerate(self.matrix[vertex_idx, :]):             if r == 1 and i in group:                 ones_in_row += 1          for i, r in enumerate(self.matrix[:, vertex_idx]):             if r == 1 and i in group:                 ones_in_column += 1          return max(ones_in_row, ones_in_column)      def _strong_connected_components(self):         return tarjan(matr2dct(self.matrix))      def initialize(self):         self.scc = self._strong_connected_components()         self.graph = nx.DiGraph(self.matrix)      def get_cycles(self):         result = []         nodes = set()         lst = list(nx.simple_cycles(self.graph))         lst.sort(key=len, reverse=True)         for cycle in lst:             if all((n not in nodes for n in cycle)):                 result.append(cycle)                 for n in cycle:                     nodes.add(n)         return result      def get_component_by_vertex(self, vertex_idx):         cycles = self.get_cycles()         for cycle in cycles:             if vertex_idx in cycle:                 return cycle      def find_to_delete(self, cycle):         dct = dict(zip(cycle, [self.cycling_order(c) for c in cycle]))         inv_dct = {}         for k, v in dct.items():             inv_dct[v] = [key for key in dct if dct[key] == v]         maxKey = max(inv_dct.keys())         candidates = inv_dct[maxKey]         if len(candidates) == 1:             return candidates[0], False         group = self.get_component_by_vertex(candidates[0])         return min(candidates, key=group.index), True      def break_loops(self):         N = self.matrix.shape[0]         self.translator = dict(zip(list(range(N)), list(range(N))))         for cycle in self.get_cycles():             num, reason = self.find_to_delete(cycle)             self.equations_indices.append(num)             print(f'Для цикла {cycle} будет удалена вершина: {num}. Причина: Максимальный цикловый порядок' + ', первым встретился в алгоритме Тарьяна' * reason)             self.matrix = np.delete(self.matrix, num, axis=0)             self.matrix = np.delete(self.matrix, num, axis=1)             del self.translator[num]             for k in self.translator:                 if k > num:                     self.translator[k] -= 1         self.inverted_translator = {v: k for k, v in self.translator.items()}         self.initialize()      def get_answer(self):         print(f'Индексы уравнений: {self.equations_indices}')         for sublist in self.scc:             if len(sublist) > 1:                 print(f'{sublist} не подстановка, а видимо уравнение. Добавить в eq_indices?')                 continue             self.substitution_indices.append(self.inverted_translator[sublist[0]])         print(f'Индексы(!) подстановок: {self.substitution_indices}')         print(f'Индексы(!) подстановок, которые действительно есть в системе: {[i for i in self.substitution_indices if i in self.in_system]}')

Дополнительный пример

Под конец хотелось бы обратить внимание на другой способ инициализации системы:

  tsk = SubstitutionDetector(8)   tsk.add_equation(1, [2, 5])   tsk.add_equation(3, [1, 6])   tsk.add_equation(2, [4, 7])   tsk.add_equation(4, [3, 8])   tsk.add_equation(6, [7])    tsk.initialize()   tsk.break_loops()   tsk.get_answer()

Для цикла [0, 1, 3, 2] будет удалена вершина: 0. Причина: Максимальный цикловый порядок, первым встретился в алгоритме Тарьяна.

Индексы уравнений: [0] Индексы(!) подстановок: [6, 5, 2, 7, 3, 1, 4] Индексы(!) подстановок, которые действительно есть в системе: [5, 2, 3, 1]

В этом случае система не содержит уравнений, описывающие переменные x_7, x_8, x_5 , но алгоритм работает, упорядочивая и классифицируя существующие уравнения.

ссылка на оригинал статьи https://habr.com/ru/articles/831950/