在python中使用有限差分的雅可比矩阵的行列式

Question

我正在尝试计算雅可比矩阵的行列式，并根据函数 x 和 y 评估该行列式何时为零。我正在尝试实现有限差分的简单方法，但结果似乎不正确。当然有更好的方法可以做到这一点，但我已经搜索了几天，但我没有得到解决方案

PS：雅可比矩阵应该是这样的

还有我正在使用的代码

import numpy as np

def x(theta1, theta2, w, h, L1, L2):
    sint1 = np.sin(theta1)
    cost1 = np.cos(theta1)
    sint2 = np.sin(theta2)
    cost2 = np.cos(theta2)

    i1 = L1 * (cost1 + cost2) + w
    j1 = L1 * (sint1 - sint2) - h
    D = np.sqrt((L1*(cost2-cost1)+w)**2+(L1*(sint2-sint1)+h)**2)
    a = (0.25)*np.sqrt((4*L2**2-D**2)*D**2)

    return i1/2 + 2*j1*a/(D**2)

def y(theta1, theta2, w, h, L1, L2):
    sint1 = np.sin(theta1)
    cost1 = np.cos(theta1)
    sint2 = np.sin(theta2)
    cost2 = np.cos(theta2)

    i2 = L1 * (sint1 + sint2) + h
    j2 = L1 * (cost1 - cost2) - w
    D = np.sqrt((L1*(cost2-cost1)+w)**2+(L1*(sint2-sint1)+h)**2)
    a = (0.25)*np.sqrt((4*L2**2-D**2)*D**2)

    return i2/2 - 2*j2*a/(D**2)


def     det_jacobiano(theta1, theta2, eps, w, h, L1, L2):    

    dxdt1 = (x(theta1+e    ps, theta2, w, h, L1, L2)-x(theta1, theta2, w,     h, L1, L2)) / (eps)    
    dxdt2 = (x(theta1,     theta2+eps, w, h, L1, L2)-x(theta1, theta2, w,     h, L1, L2)) / (eps)    
    dydt1 = (y(theta1+e    ps, theta2, w, h, L1, L2)-y(theta1, theta2, w,     h, L1, L2)) / (eps)    
    dydt2 = (y(theta1, theta2+eps, w, h, L1, L2)-y(theta1, theta2, w,     h, L1, L2)) / (eps)

    return dxdt1*dydt2-dxdt2*dydt1


def singular    idades(theta1_min,theta1_max, theta2_min,theta2_max, n,     eps,tol, w,     h, L1, L2):    
    x_s = []    
    y_s = []    
    theta1 = np.linspace(theta1_    min,theta1_max,n)    
    theta2 = np.linspace(theta2_    min,theta2_max,n)    
    for i in range(len(theta1)):    
        for j in range(len(theta2)):    
            det_jac = det_jacobiano(theta1[i], theta2[j], eps, w, h,     L1, L2)
            if det_jac<tol and det_jac>-tol:
                x_s.append(x(theta1[i], theta2[j], w, h, L1, L2))
                y_s.append(y(theta1[i], theta2[j], w, h, L1, L2))
    return x_s, y_s

x_s, y_s = singularidades(0,2*np.pi,0,2*np.pi,100,1e-8,0.0001, 0, 0, 100, 100)

Answer 1

Numdifftools 应该这样做。它还提供其他 difftools。