Ctrl+K

Нелинейные системы уравнений

Нелинейные системы уравнений#

Пусть имеем систему нелинейных уравнений

\[\begin{split} \left\{\begin{array}{l} f_1\left(x_1, x_2, \ldots, x_n\right)=0 \\ f_2\left(x_1, x_2, \ldots, x_n\right)=0 \\ ...\\ f_n\left(x_1, x_2, \ldots, x_n\right)=0 \end{array} \quad \Leftrightarrow \quad \vec{F}(x)=0\right. \end{split}\]

Хотим найти все её корни.

Сведение к поиску минимума#

Перепишем её в виде одной функции \(n\) переменных:

\[ \Psi(x)=\sum_{i=1}^n\left(f_n\right)^2=0 \]

А тогда решение системы эквивалентно поиску глобальных минимумов этой функции:

\[ \min \Psi(x) = \, ? \]

А дальше глобальная минимизация ;) - так что такой метод не особо применим в сложных задачах.

Метод простой итерации (релаксации)#

Аналогично параграфу про МПИ для функции одной переменной:

\[\begin{split} \begin{aligned} &\vec{F}(x)=0 \iff \vec{\Phi}(x)=x \\ \\ &x_{n+1}=\vec{\Phi}\left(x_n\right) \end{aligned} \end{split}\]

Заметим, что домножение на невырожденную матрицу \(\Lambda\) (которая также может зависеть от \(x\)) не влияет на корни. Таким образом, получаем обобщение метода релаксации для нескольких переменных:

\[ \vec{F}(x)=0 \iff \Lambda \vec{F}(x)=0 \]

и тогда \(\vec{\Phi}(x)\) можно выбрать как

\[ \vec{\Phi}(x)=x - \Lambda \vec{F}(x) \]

и решать прежний итерационный процесс:

Шаг метода релаксации

\[ x_{n+1}=\vec{\Phi}\left(x_n\right) \]

где \(\vec{\Phi}(x)=x - \Lambda(x) \vec{F}(x)\) с произвольной \(\Lambda(x)\) такой что \(\operatorname{det} \Lambda(x) \neq 0\)

Достаточное условие сходимости МПИ

Достаточным условием сжимаемости отображения будет:

\[ \max _{x \in G} \max _i \sum_{j=1}^n\left|\frac{\partial \varphi_i(x)}{\partial x_j}\right| \leq q<1 \]

или, что эквивалентно, какая-либо норма матрицы Якоби \(\vec{\Phi}(x)\) должна быть меньше единицы. \(\Lambda(x)\) надо выбирать именно из этого условия.

Метод Зейделя#

\[ \vec{F}(x)=0 \iff \vec{\Phi}(x)=x \]

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

Шаг Зейделя

\[\begin{split} \begin{aligned} &x_1^{n+1}=\varphi_1\left(x_1^n, x_2^n, ..., x_n^n\right) \\ \\ &x_2^{n+1}=\varphi_2\left(x_1^{n+1}, x_2^n, \ldots, x_n^n\right) \\ \\ &x_3^{n+1}=\varphi_3\left(x_1^{n+1}, x_2^{n+1}, x_3^{n} \ldots, x_n^n\right) \\ \\ &... \end{aligned} \end{split}\]

Ожидается улучшение сходимости по сравнению с обычным МПИ.

Метод Ньютона#

Является частным случаем метода релаксации для системы при выборе \(\Lambda (x) = J^{(-1)}(x)\) - обратная матрицы Якоби вектора уравнений \(\vec{F}\). Тогда формула для итерации будет:

\[ x_{n+1}=x_n-J^{(-1)}\left(x_n\right) \cdot \vec{F}\left(x_n\right) \]

Однако вычислять эту обратную матрицу Якоби вовсе необязательно. Заметим, что эта формула эквивалентна следующей системе линейных уравнений:

Шаг метода Ньютона для системы уравнений

\[ J\left(x_n\right) \Delta x_{n+1}=-\vec{F}\left(x_n\right) \]

Решая эту СЛАУ относительно вектора \(\Delta x_{n+1}\), мы и получаем прибавку на текущей итерации.

Пример. Проиллюстрируем этот метод на решении следующей системы:

\[\begin{split} \left\{\begin{array}{l} x_1^2-x_2+x_1 \cdot \cos \left(\pi x_1\right)=0 \\ x_1 \cdot x_2+e^{-x_2}-\frac{1}{x_1}=0 \end{array}\right. \end{split}\]

Источник.

Для вычисления матрицы Якоби воспользуемся sympy.

# Не хотим вычислять матрицу Якоби ручками - символьные вычисления спешат на помощь
import sympy as smp

# Определение переменных
x1, x2 = smp.symbols('x1 x2')

# Определение уравнений
eq1 = x1**2 - x2 + x1*smp.cos(smp.pi*x1)
eq2 = x1*x2 + smp.exp(-x2) - 1/x1

# Создание списка уравнений
equations = [eq1, eq2]

# Нахождение Якобиана
Jacobian_matrix = smp.Matrix([[smp.diff(eq, var) for var in (x1, x2)] for eq in equations])

# Вывод уравнений и матрицы Якоби
print("Система уравнений:")
for eq in equations:
    display(eq)

print("\nМатрица Якоби:")
display(Jacobian_matrix)

Система уравнений:
\[\displaystyle x_{1}^{2} + x_{1} \cos{\left(\pi x_{1} \right)} - x_{2}\]
\[\displaystyle x_{1} x_{2} + e^{- x_{2}} - \frac{1}{x_{1}}\]
Матрица Якоби:
\[\begin{split}\displaystyle \left[\begin{matrix}- \pi x_{1} \sin{\left(\pi x_{1} \right)} + 2 x_{1} + \cos{\left(\pi x_{1} \right)} & -1\\x_{2} + \frac{1}{x_{1}^{2}} & x_{1} - e^{- x_{2}}\end{matrix}\right]\end{split}\]
import numpy as np

def F(x): # Так то и это можно сразу же из sympy сделать, но для наглядности просто перепишем
    return np.array(
        [x[0]**2 - x[1] + x[0]*np.cos(np.pi*x[0]),
         x[0]*x[1] + np.exp(-x[1]) - x[0]**(-1)])

def J(x): # Так то и это можно сразу же из sympy сделать, но для наглядности просто перепишем
    return np.array(
        [[2*x[0] + np.cos(np.pi*x[0]) - np.pi*x[0]*np.sin(np.pi*x[0]), -1],
         [x[1] + x[0]**(-2), x[0] - np.exp(-x[1])]])

def Newton_system(F, J, x, eps):
    """
    Solve nonlinear system F=0 by Newton's method.
    J is the Jacobian of F. Both F and J must be functions of x.
    At input, x holds the start value. The iteration continues
    until ||F|| < eps.
    """
    F_value = F(x)
    F_norm = np.linalg.norm(F_value, ord=2)  # l2 norm of vector
    iteration_counter = 0
    while abs(F_norm) > eps and iteration_counter < 100:
        delta = np.linalg.solve(J(x), -F_value)
        x = x + delta
        F_value = F(x)
        F_norm = np.linalg.norm(F_value, ord=2)
        iteration_counter += 1

    # Here, either a solution is found, or too many iterations
    if abs(F_norm) > eps:
        iteration_counter = -1
    return x, iteration_counter

expected = np.array([1, 0])
tol = 1e-4
x, n = Newton_system(F, J, x=np.array([2., -1.]), eps=0.0001)
print(x, n)
error_norm = np.linalg.norm(expected - x, ord=2)
assert error_norm < tol, 'norm of error =%g' % error_norm
print('norm of error =%g' % error_norm)

[ 1.00000006e+00 -1.00943962e-06] 4
norm of error =1.01115e-06
Содержание