用Numdifftools计算标量场的Hessian矩阵；与MATLAB相比，结果出人意料

Question

我试图用“Numdifftools”计算标量场的Hessian矩阵。然而，我在Hessian矩阵的一个元素中得到了一些意想不到的结果

这是我的标量场。它是两个变量（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

以下是我的hessian计算例程：

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])

我知道我的Hessian矩阵，对于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