牛顿-拉夫逊法求解闭矢量环二元方程组无法收敛到0.001精度的问题求助

我写了一个用牛顿-拉夫逊法求解二元方程组根的程序，目的是求解闭矢量环的两个未知量。但问题是，只有当误差容限设为0.01及更大时程序才能正常工作，当容限设为比如0.001时，程序会陷入振荡，似乎永远无法收敛。我觉得这种行为很奇怪，毕竟它能达到0.01的精度，说明初始猜测应该已经进入了收敛半径，按道理后续的近似不应该有问题...

以下是我的代码：

import numpy as np
import random as rnd

TOLERANCE = 0.01
ITER_LIMIT = 50

track_x = 0
track_y = 0
track_incline = np.deg2rad(60)

r41 = 12.2/2
r1 = 28.88 - 10.95
r2 = 22.65 - track_x
r26 = 8/2 - track_y
r3 = 8
r4 = 12.2

theta41 = np.pi
theta3 = np.deg2rad(30)

def f1(x, y):
    return r41*np.sin(theta41) + r1*np.cos(theta41) + r2*np.sin(y) + r26*np.cos(y) - x*np.sin(y + track_incline) + r3*np.sin(theta3)

def f2(x, y):
    return r41*np.cos(theta41) - r1*np.sin(theta41) + r2*np.cos(y) - r26*np.sin(y) - x*np.cos(y + track_incline) + r3*np.cos(theta3) + r4

def jacobian_mat(x, y):
    return np.array([
        [-np.sin(y + track_incline), r2*np.cos(y) - r26*np.sin(y) - x*np.cos(y + track_incline)],
        [-np.cos(y + track_incline), -r2*np.sin(y) - r26*np.cos(y) + x*np.sin(y + track_incline)]
    ])

def compute_newton_raphson(x0, y0):
    for i in range(ITER_LIMIT):
        result = -np.matmul(np.array([f1(x0, y0), f2(x0, y0)]), np.linalg.inv(jacobian_mat(x0, y0))) + np.array([x0, y0])
        x = result[0]
        y = result[1]

        print(x, y)

        if max(np.abs(x - x0), np.abs(y - y0)) <= TOLERANCE:
            return (x, y)
        else:
            x0 = x
            y0 = y

    return None

while True:
    x0 = rnd.randrange(-100, 100)
    y0 = rnd.randrange(-100, 100)

    res = compute_newton_raphson(x0, y0)
    if res != None:
        print(f"RESULT:\nX = {res[0]}mm\nY = {np.rad2deg(res[1])}deg")
        break

我准备了问题的建模草图和角度标注图，或许能帮大家理解这个问题的背景。

我尝试过调整迭代次数，但完全没用；也仔细检查了f1和f2的构造是否正确，确认没有问题。

希望能得到大家的任何见解，非常感谢！

备注：内容来源于stack exchange，提问作者Rapolas Daugelis