Не могу разобраться с методами наискорейшего спуска

Question

Пишу программу, которая вычисляет корни нелинейных уравнений при помощи метода наискорейшего спуска. У меня проблемы также возникли с поисками производной (приложу скриншот, чтобы было понятнее, о чём идёт речь). Хотелось бы использовать методы numpy и sumpy, но так как времени уже не так много до экзамена остаётся, поэтому написал код на скорую руку.

Первый случай x1 = -0.9, x2 = 0.7 (корни должны быть с точностью 10^-6: x1 = -0,947718 и x2 = 0,713550, но вместо этого совсем другой корень находит).

Второй случай x1 = -0.1, x2 = 2.2 (корни должны быть с точностью 10^-6: x1 = -0,110086 и x2 = 2,222477, как это было в методе Ньютона).

Система нелинейных уравнений:

Графическая локализация корня уравнения:

Теорию можно увидеть тут:

Метод наискорейшего спуска, страница 1

Метод наискорейшего спуска, страница 2

Мои попытки решить это уравнение:

Код программы, начальная точка (-0.9, 0.7):

import math

# первый корень
x1k = -0.9
x2k = 0.7

# второй корень
##x1k = -0.1
##x2k = 2.2

alpha = 0.25

##    fx1_x1k_x2k = round(2 * (math.cos(x1k + 1.5) * (math.sin(x1k + 1.5) + 2.9 - x2k) + x1k + math.cos(x2k - 2)), 6)
##    fx2_x1k_x2k = round(-2 * ((math.cos(x2k - 2) * math.sin(x2k - 2)) - x2k + math.sin(x1k + 1.5) + 2.9 + (x1k * math.sin(x1k + 1.5))), 6)

for k in range(1, 39):
    fx1_x1k_x2k = round((2 * math.sin(x1k + x2k) * math.cos(x1k + x2k)) + (3.38 * x1k) - (2.6 * math.sin(x1k + x2k)) - (2.6 * x1k * math.cos(x1k + x2k)) - (2 * math.cos(x1k + x2k)) + 2.6, 6)
    fx2_x1k_x2k = round(4 * x1k**3 + 0.16 * x2k**3 + (0.8 * x2k * x1k**2) + (0.4 * x2k**2 * 2 * x1k) - 4 * x1k - 0.8 * x2k, 6)
    
    x1_current = round(x1k - alpha * fx1_x1k_x2k, 6)
    x2_current = round(x2k - alpha * fx2_x1k_x2k, 6)

    print(f'------------------------------------ k = {k} ------------------------------------')
    
    print(f'x1_current = {x1_current}')
    print(f'x2_current = {x2_current}')

    print(f'|x1 - x1_current| = {round(abs(x1k - x1_current), 6)}')

    x1k = x1_current
    x2k = x2_current

Результат нахождения первого корня, вместо (-0.947718, 0.713550) получаем (-0.110084, 2.222471) при k = 38:

------------------------------------ k = 1 ------------------------------------
x1_current = -0.904586
x2_current = 0.63008
|x1 - x1_current| = 0.004586
------------------------------------ k = 2 ------------------------------------
x1_current = -0.920633
x2_current = 0.550413
|x1 - x1_current| = 0.016047
------------------------------------ k = 3 ------------------------------------
x1_current = -0.950972
x2_current = 0.475969
|x1 - x1_current| = 0.030339
------------------------------------ k = 4 ------------------------------------
x1_current = -0.996371
x2_current = 0.432886
|x1 - x1_current| = 0.045399

...

------------------------------------ k = 38 ------------------------------------
x1_current = -0.110084
x2_current = 2.222471
|x1 - x1_current| = 1e-06

Код программы, начальная точка (-0.1, 2.2):

import math

# первый корень
##x1k = -0.9
##x2k = 0.7

# второй корень
x1k = -0.1
x2k = 2.2

alpha = 0.25

##    fx1_x1k_x2k = round(2 * (math.cos(x1k + 1.5) * (math.sin(x1k + 1.5) + 2.9 - x2k) + x1k + math.cos(x2k - 2)), 6)
##    fx2_x1k_x2k = round(-2 * ((math.cos(x2k - 2) * math.sin(x2k - 2)) - x2k + math.sin(x1k + 1.5) + 2.9 + (x1k * math.sin(x1k + 1.5))), 6)

for k in range(1, 39):
    fx1_x1k_x2k = round((2 * math.sin(x1k + x2k) * math.cos(x1k + x2k)) + (3.38 * x1k) - (2.6 * math.sin(x1k + x2k)) - (2.6 * x1k * math.cos(x1k + x2k)) - (2 * math.cos(x1k + x2k)) + 2.6, 6)
    fx2_x1k_x2k = round(4 * x1k**3 + 0.16 * x2k**3 + (0.8 * x2k * x1k**2) + (0.4 * x2k**2 * 2 * x1k) - 4 * x1k - 0.8 * x2k, 6)
    
    x1_current = round(x1k - alpha * fx1_x1k_x2k, 6)
    x2_current = round(x2k - alpha * fx2_x1k_x2k, 6)

    print(f'------------------------------------ k = {k} ------------------------------------')
    
    print(f'x1_current = {x1_current}')
    print(f'x2_current = {x2_current}')

    print(f'|x1 - x1_current| = {round(abs(x1k - x1_current), 6)}')

    x1k = x1_current
    x2k = x2_current

Результат нахождения второго корня, получаем (-0.110084, 2.222471) при k = 21:

------------------------------------ k = 1 ------------------------------------
x1_current = -0.106128
x2_current = 2.20748
|x1 - x1_current| = 0.006128
------------------------------------ k = 2 ------------------------------------
x1_current = -0.105683
x2_current = 2.212223
|x1 - x1_current| = 0.000445
------------------------------------ k = 3 ------------------------------------
x1_current = -0.108148
x2_current = 2.215606
|x1 - x1_current| = 0.002465
------------------------------------ k = 4 ------------------------------------
x1_current = -0.108136
x2_current = 2.217791
|x1 - x1_current| = 1.2e-05

...

------------------------------------ k = 21 ------------------------------------
x1_current = -0.110084
x2_current = 2.222471
|x1 - x1_current| = 1e-06

Answer 1

Вот код, который использует sympy для поиска производной, и numpy для простоты кода:

import numpy as np
import sympy as sp


def gradient_descent_system(equations, variables, start_point, learning_rate=0.01, tolerance=1e-6, max_iter=1000):
    # Функция которую нужно минимизировать
    F = sum(eq ** 2 for eq in equations)

    # Производные функции
    gradient = [sp.diff(F, var) for var in variables]
    gradient_func = sp.lambdify(variables, gradient, 'numpy')

    x = start_point.astype(float)

    for i in range(max_iter):
        # Находим градиент
        grad = np.array(gradient_func(*x))

        # Двигаем точку
        x_new = x - learning_rate * grad

        # Если сдвинулись мало значит нашли корень
        if np.linalg.norm(x_new - x) < tolerance:
            break
        x = x_new

    return x


if __name__ == "__main__":
    x1, x2 = sp.symbols('x1 x2')

    eq1 = sp.sin(x1 + x2) - 1.3 * x1 - 1
    eq2 = x1 * x1 + 0.2 * x2 * x2 - 1

    start_point = np.array([-0.9, 0.7])
    solution = gradient_descent_system([eq1, eq2], [x1, x2], start_point, learning_rate=0.1)
    print(f"Решение системы из {start_point}: {solution}")

    start_point = np.array([-0.1, 2.2])
    solution = gradient_descent_system([eq1, eq2], [x1, x2], start_point, learning_rate=0.1)
    print(f"Решение системы из {start_point}: {solution}")

Результат:

Решение системы из [-0.9  0.7]: [-0.9477161   0.71355646]
Решение системы из [-0.1  2.2]: [-0.11008479  2.222472  ]

БЛОГ НА HUSL

Не могу разобраться с методами наискорейшего спуска

Ответы (1 шт):