Гессенская матрица вычисления скалярного поля с использованием Numdifftools; неожиданные результаты по сравнению с MATLAB

Question

Я пытаюсь вычислить гессенскую матрицу скалярного поля, используя "Numdifftools". Я, однако, получаю некоторые неожиданные результаты в одном из элементов моей гессианской матрицы.

Вот мое скалярное поле. Это функция двух переменных (x и y), заданных x [0] и x [1] соответственно, которая возвращает lambda_2.

def lambda2Field(x):
    out = solve_ivp(doubleGyreVar, t_span=(0, T), y0=[x[0], x[1], 1, 0, 0, 1],t_eval=[0, T], rtol = 1e-6, atol=1e-6)
    output = out.y
    J = output[2:,-1].reshape(2,2,order="F")
    CG = np.matmul(J.T , J)
    lambdas, xis = np.linalg.eig(CG)
    xi_1 = xis[:,np.argmin(np.abs(lambdas))] 
    xi_2 = xis[:,np.argmax(np.abs(lambdas))]
    lambda_1 = np.min(np.abs(lambdas))
    lambda_2 = np.max(np.abs(lambdas))
    return lambda_2

Вот моя рутина гессенских вычислений:

def pointinU0(xval, yval):
    out = solve_ivp(doubleGyreVar, t_span=(0, T), y0=[xval, yval, 1, 0, 0, 1],t_eval=[0, T/2, T], rtol = 1e-10, atol=1e-10)
    output = out.y
    J = output[2:,-1].reshape(2,2,order="F")
    CG = np.matmul(J.T , J)
    lambdas, xis = np.linalg.eig(CG)
    xi_1 = xis[:,np.argmin(np.abs(lambdas))] 
    xi_2 = xis[:,np.argmax(np.abs(lambdas))]
    lambda_1 = np.min(np.abs(lambdas))
    lambda_2 = np.max(np.abs(lambdas))
    hessian = nd.Hessian(lambda2Field)([xval,yval])

Я знаю, что моя матрица Гессе должна при x = [0,0] быть равна

1.0e+20 *

  [-0.012176339633369  -0.000000000000000
  -0.000000000000000  -2.375233755791394]

Как вычислено с MATLAB. Тем не менее, мой код дает следующее:

array([[-1.21751109e+18,  0.00000000e+00],
       [ 0.00000000e+00, -9.79794699e+16]])

Интересно, что отличается только нижний правый компонент, и я понятия не имею, откуда это может появиться. Я попытался повторить следующий код MATLAB для определения поля lambda_2.

function lam2=lambda(x);
        T=[0 20];
        inneropt = odeset('AbsTol',1e-6, 'RelTol', 1e-6);

        id = eye(length(x));
        [~,xrr] = ode113(@variational3,[T(1) T(2)],[x ; id(:)],inneropt);
        JAC=reshape(xrr(end, 3:6), 2, 2);
        CG=JAC'*JAC;
        [E,L]=eig(CG);


            if abs(L(1,1)) < abs(L(2,2))
                l1=abs(L(1,1));
                l2=abs(L(2,2));
            else
                l1=abs(L(2,2));
                l2=abs(L(1,1));
            end
            lam2=l2;
end