Вычисление корней уравнений методом Ньютона с помощью якобиана системы

Пытаюсь написать программу, которая вычисляет корни уравнений методом Ньютона (с помощью нахождения якобиана системы).

Текст задания с решением (файл обновлён, но пока не всё решил до конца):

  1. Локализуйте корни системы уравнений графически.

sin(x1 + x2) - 1,3x1 - 1 = 0

(x1)^2 + (0,2x2)^2 - 1 = 0

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

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

  1. Найдите с точностью ε = 10^-6 все корни системы нелинейных уравнений, используя методы Ньютона и наискорейшего спуска.

Формулы, которые я использую для вычислений:

Формулы для вычислений

Программа на Python:

import math

# первый корень
x1 = -0.9
x2 = 0.7

for k in range(10):
    print(f'------------------------------------ k = {k + 1} ------------------------------------')
    
    a = round(math.cos(x1 + x2), 6)
    b = 1.3
    c = round(2 * x1, 6)
    d = round(0.4 * x2, 6)

    print(f'x1 = {x1} (при k = {k})')
    print(f'x2 = {x2} (при k = {k})')

    print(f'x1 + x2 = {round(x1 + x2, 6)}')
    
    A = [[a, b], [c, d]]
    print(f'A = {A}')
    
    determinant = round(A[0][0] * A[1][1] - A[0][1] * A[1][0], 6)
    print(f'det A = {determinant}')
    
    A_union = [[((-1)**(1 + 1)) * A[1][1], ((-1)**(1 + 2)) * A[1][0]],
             [((-1)**(2 + 1)) * A[0][1], ((-1)**(2 + 2)) * A[0][0]]]
    print(f'A_union = {A_union}')

    A_union_T = [[0, 0], [0, 0]]
    for j in range(len(A)):
        for i in range(len(A)):
            A_union_T[i][j] = A_union[j][i]
    print(f'A_union_T = {A_union_T}')

    A_reversed = [[round((1/determinant) * A_union_T[0][0], 6), round((1/determinant) * A_union_T[0][1], 6)],
                  [round((1/determinant) * A_union_T[1][0], 6), round((1/determinant) * A_union_T[1][1], 6)]]
    print(f'A_reversed = {A_reversed}')

    fx0 = [round(math.sin(x1 + x2) - (1.3 * x1) - 1, 6), round(x1**2 + (0.2 * x2**2) - 1, 6)]
    print(f'fx0 = {fx0}')

    print(f'A_reversed[0][0] * fx0[0] + A_reversed[0][1] * fx0[1] {round(A_reversed[0][0] * fx0[0] + A_reversed[0][1] * fx0[1], 6)}')
    print(f'A_reversed[1][0] * fx0[0] + A_reversed[1][1] * fx0[1] {round(A_reversed[1][0] * fx0[0] + A_reversed[1][1] * fx0[1], 6)}')

    x1_current = round(x1 - (A_reversed[0][0] * fx0[0] + A_reversed[0][1] * fx0[1]), 6)
    x2_current = round(x2 - (A_reversed[1][0] * fx0[0] + A_reversed[1][1] * fx0[1]), 6)
    print(f'x1_current = {x1_current}')
    print(f'x2_current = {x2_current}')

    print(f'||x(k) - x(k + 1)|| = {abs(round(x1 - x1_current, 6))}')

    x1 = x1_current
    x2 = x2_current

UPD: В первом случае, когда начальное приближение x1 = -0.9 и x2 = 0.7, то корень мне удалось найти при k = 10 (x1 = -0.947577, x2 = 0.714488). Но во втором случае когда я пытаюсь корень найти при начальном приближении x1 = -0.1 и x2 = 2.2, то не получается такой же программой найти решение. По идее, должно быть x1 = -0.1 с чем-то там и x2 = 2.2 с чем-то там.

------------------------------------ k = 1 ------------------------------------
x1 = -0.1 (при k = 0)
x2 = 2.2 (при k = 0)
x1 + x2 = 2.1
A = [[-0.504846, 1.3], [-0.2, 0.88]]
det A = -0.184264
A_union = [[0.88, 0.2], [-1.3, -0.504846]]
A_union_T = [[0.88, -1.3], [0.2, -0.504846]]
A_reversed = [[-4.775757, 7.055095], [-1.085399, 2.739797]]
fx0 = [-0.006791, -0.022]
A_reversed[0][0] * fx0[0] + A_reversed[0][1] * fx0[1] -0.12278
A_reversed[1][0] * fx0[0] + A_reversed[1][1] * fx0[1] -0.052905
x1_current = 0.02278
x2_current = 2.252905
||x(k) - x(k + 1)|| = 0.12278
------------------------------------ k = 2 ------------------------------------
x1 = 0.02278 (при k = 1)
x2 = 2.252905 (при k = 1)
x1 + x2 = 2.275685
A = [[-0.647949, 1.3], [0.04556, 0.901162]]
det A = -0.643135
A_union = [[0.901162, -0.04556], [-1.3, -0.647949]]
A_union_T = [[0.901162, -1.3], [-0.04556, -0.647949]]
A_reversed = [[-1.401202, 2.021349], [0.07084, 1.007485]]
fx0 = [-0.26793, 0.015635]
A_reversed[0][0] * fx0[0] + A_reversed[0][1] * fx0[1] 0.407028
A_reversed[1][0] * fx0[0] + A_reversed[1][1] * fx0[1] -0.003228
x1_current = -0.384248
x2_current = 2.256133
||x(k) - x(k + 1)|| = 0.407028

...

------------------------------------ k = 10 ------------------------------------
x1 = -0.947328 (при k = 9)
x2 = 0.71615 (при k = 9)
x1 + x2 = -0.231178
A = [[0.973397, 1.3], [-1.894656, 0.28646]]
det A = 2.741892
A_union = [[0.28646, 1.894656], [-1.3, 0.973397]]
A_union_T = [[0.28646, -1.3], [1.894656, 0.973397]]
A_reversed = [[0.104475, -0.474125], [0.691003, 0.355009]]
fx0 = [0.002402, 5e-06]
A_reversed[0][0] * fx0[0] + A_reversed[0][1] * fx0[1] 0.000249
A_reversed[1][0] * fx0[0] + A_reversed[1][1] * fx0[1] 0.001662
x1_current = -0.947577
x2_current = 0.714488
||x(k) - x(k + 1)|| = 0.000249

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