计算多变量函数在Python中的区间之间取值的有效方法

Question

我有两个确定机器人逆运动学的多元方程。这些方程取决于变量 theta1 和 theta2（其他变量是几何常数）

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_jacobian(theta1, theta2, w, h, L1, L2,eps):
    dxdt1 = (-x(theta1+eps, theta2, w, h, L1, L2)+4*x(theta1, theta2, w, h, L1, L2)-3*x(theta1-eps, theta2, w, h, L1, L2))/(2*eps)
    dxdt2 = (-x(theta1, theta2+eps, w, h, L1, L2)+4*x(theta1, theta2, w, h, L1, L2)-3*x(theta1, theta2-eps, w, h, L1, L2))/(2*eps)
    dydt1 = (-y(theta1+eps, theta2, w, h, L1, L2)+4*y(theta1, theta2, w, h, L1, L2)-3*y(theta1-eps, theta2, w, h, L1, L2))/(2*eps)
    dydt2 = (-y(theta1, theta2+eps, w, h, L1, L2)+4*y(theta1, theta2, w, h, L1, L2)-3*y(theta1, theta2-eps, w, h, L1, L2))/(2*eps)  
    return dxdt1,dxdt2,dydt1,dydt2

针对属于某个区间的 theta1 和 theta2 的值进行评估

theta1 = np.linspace(theta1_min,theta1_max,n)
theta2 = np.linspace(theta2_min,theta2_max,n)
theta1, theta2 = np.meshgrid(theta1,theta2)

我想知道是否有一种有效的方法（使用 numpy 数组）来计算 x 和 y 的值，其中行列式的值介于 -tol 和 tol 之间（tol=1e-08）。目前我正在使用两个嵌套的for循环，但速度很慢

我用for循环写了一个函数，但是很慢

def singularidades(theta1_min,theta1_max, theta2_min,theta2_max, n,tol, w, h, L1, L2,eps):
    x_s = []
    y_s = []
    theta1_s = []
    theta2_s = []
    det = []
    theta1 = np.linspace(theta1_min,theta1_max,n)
    theta2 = np.linspace(theta2_min,theta2_max,n)
    theta1, theta2 = np.meshgrid(theta1,theta2)
    det_jac = det_jacobiano(theta1,theta2,w,h,L1,L2,eps)
    for i in range(n):
        for j in range(n):
            if (g_tol[i,j] and l_tol[i,j]):
                x_s.append(x(theta1[i,j], theta2[i,j], w, h, L1, L2))
                y_s.append(y(theta1[i,j], theta2[i,j], w, h, L1, L2))
                theta1_s.append(theta1[i,j])
                theta2_s.append(theta2[i,j])
                det.append(det_jac[i,j])
   return x_s,y_s,theta1_s,theta2_s,det,(g_tol and l_tol)

编辑：我已修改 det_jacobian 函数以将其与 scipy.optimize.root 一起使用

def det_jacobiano(theta, w, h, L1, L2,eps):
    theta1,theta2 = theta
    dxdt1 = (-x(theta1+eps, theta2, w, h, L1, L2)+4*x(theta1, theta2, w, h, L1, L2)-3*x(theta1-eps, theta2, w, h, L1, L2))/(2*eps)
    dxdt2 = (-x(theta1, theta2+eps, w, h, L1, L2)+4*x(theta1, theta2, w, h, L1, L2)-3*x(theta1, theta2-eps, w, h, L1, L2))/(2*eps)
    dydt1 = (-y(theta1+eps, theta2, w, h, L1, L2)+4*y(theta1, theta2, w, h, L1, L2)-3*y(theta1-eps, theta2, w, h, L1, L2))/(2*eps)
    dydt2 = (-y(theta1, theta2+eps, w, h, L1, L2)+4*y(theta1, theta2, w, h, L1, L2)-3*y(theta1, theta2-eps, w, h, L1, L2))/(2*eps)  
    return dxdt1*dydt2 - dxdt2*dydt1

我正在尝试使用

initial_guess = [2.693, 0.4538]
result = optimize.root(det_jacobiano, initial_guess,tol=1e-8,args=(20,0,100,100,1e-10),method='lm')

但我得到了错误：

TypeError: Improper input: N=2 must not exceed M=1

Answer 1

您不需要循环。您的函数可以使用 numpy 数组以及单个值：

def f(x,y):
    return np.sin(x + y) / np.sqrt(x**2 + y**2)

x = [0.1, 0.2, 0.3, 0.4, 0.5]
y = [0.1, 0.2, 0.3, 0.4, 0.5]

print(f(x, y))

将返回：

[1.40480431, 1.37680175, 1.33087507, 1.26811839, 1.19001968]

它是每对 x 和 y 的函数值数组

Answer 2

一种方法是定义：

def f(x,y,a):
    return np.sin(x + y) / np.sqrt(x**2 + y**2) - a

其中a 可以获取区间的限制，然后使用scipy.optimize 计算此函数的根x0 和y0。

这些根对应于函数的 x 和 y 值，函数在这两个值之间返回所需的值。