Python，四杆机构角度时间图

Question

我试图绘制四杆机构输出角度的角度与时间图（下图中的角度fi4）。使用第23页https://scholar.cu.edu.eg/?q=anis/files/week04-mdp206-position_analysis-draft.pdf中的解计算该角度

我现在正试图绘制fi_4(t)图，得到了一些奇怪的结果。该图将输入角度fi2显示为蓝色，输出角度fi4显示为红色。为什么fi2会随着时间而波动？{}不应该有某种正弦曲线吗

我是不是遗漏了什么

四连杆机构：

守则：

from __future__ import division
import math
import numpy as np
import matplotlib.pyplot as plt

# Input
#lengths of links (tube testing machine actual lengths)
a = 45.5   #mm
b = 250   #mm
c = 140   #mm
d = 244.244  #mm

# Solution for fi2 being a time function, f(time) = angle
f = 16.7/60    #/s
omega = 2 * np.pi * f   #rad/s

t = np.linspace(0, 50, 100)
y = a * np.sin(omega * t)
x = a * np.cos(omega * t)

fi2 = np.arctan(y/x)

# Solution of the vector loop equation
#https://scholar.cu.edu.eg/?q=anis/files/week04-mdp206-position_analysis-draft.pdf
K1 = d/a
K2 = d/c
K3 = (a**2 - b**2 + c**2 + d**2)/(2*a*c)
A = np.cos(fi2) - K1 - K2*np.cos(fi2) + K3
B = -2*np.sin(fi2)
C = K1 - (K2+1)*np.cos(fi2) + K3
fi4_1 = 2*np.arctan((-B+np.sqrt(B**2 - 4*A*C))/(2*A))
fi4_2 = 2*np.arctan((-B-np.sqrt(B**2 - 4*A*C))/(2*A))

# Plot the fi2 time diagram and fi4 time diagram
plt.plot(t, np.degrees(fi2), color = 'blue')
plt.plot(t, np.degrees(fi4_2), color = 'red')
plt.show()

图表：

Answer 1

linespace(0, 50, 100)太快了。将其替换为：

t = np.linspace(0, 5, 100)

其次，所有涉及裸np.arctan()的计算都是不正确的。您应该使用^{}，来确定正确的象限（与基于y/x的任何方法不同，x和y的各自符号丢失）。因此：

fi2 = np.arctan2(y, x)  # not: np.arctan(y/x)

...

fi4_1 = 2 * np.arctan2(-B + np.sqrt(B**2 - 4*A*C), 2*A)
fi4_2 = 2 * np.arctan2(-B - np.sqrt(B**2 - 4*A*C), 2*A)

在绘图上添加一些标签并显示θ_4的两种解决方案：

plt.plot(t, np.degrees(fi2) % 360, color = 'k', label=r'$θ_2$')
plt.plot(t, np.degrees(fi4_1) % 360, color = 'b', label=r'$θ_{4_1}$')
plt.plot(t, np.degrees(fi4_2) % 360, color = 'r', label=r'$θ_{4_2}$')
plt.xlabel('t [s]')
plt.ylabel('degrees')
plt.legend()
plt.show()

通过这些MOD，我们可以：

Answer 2

顺便说一句，你想看到一个令人惊讶的懒惰的方式来解决这些问题吗？比您的代码效率要低得多，但在不尝试表达解决方案的封闭形式的情况下更容易派生（例如，对于其他结构）：

from scipy.optimize import fsolve

def polar(r, theta):
    return r * np.array((np.cos(theta), np.sin(theta)))

def f(th34, th2):
    th3, th4 = th34  # solve simultaneously for theta_3 and theta_4
    pb_23 = polar(a, th2) + polar(b, th3)  # point B based on links a, b
    pb_14 = polar(d, 0) + polar(c, th4)  # point B based on links d, c
    return pb_23 - pb_14  # error: difference of the two

def solve(th2):
    th4_1 = np.array([fsolve(f, [0, -1.5], args=(th2_k,))[1] for th2_k in th2])
    th4_2 = np.array([fsolve(f, [0, 1.5], args=(th2_k,))[1] for th2_k in th2])
    return th4_1, th4_2

应用程序：

t = np.linspace(0, 5, 100)
th2 = omega * t
th4_1, th4_2 = solve(th2)

twopi = 2 * np.pi
np.allclose(th4_1 % twopi, fi4_1 % twopi)
# True

np.allclose(th4_2 % twopi, fi4_2 % twopi)
# True

根据机构的结构（例如5个连杆），您可能有两个以上的解决方案，当然还有更多的角度，因此您必须修改上述代码。但你明白了

请注意：fsolve迭代以找到合适的（足够接近的）解决方案，正如我所说的，它比封闭形式慢得多

更新（一些澄清/解释）：

• 函数f以两种不同的方式（通过R2-R3和通过R1-R4）计算点B的位置，并返回差值（作为向量）。我们求差为零

• 该函数接受两个参数：一个二维变量（th34，它是数组[th3, th4]）和一个参数th2；该参数在一次运行fsolve期间为常量

• 值[0, -1.5]和[0, 1.5]是th34（th3和th4）的初始化值（猜测）。我们调用fsolve两次以得到两个可能的解决方案

• 所有角度都参考你的图形。我使用th表示θ（θ，而不是φ），但我保留了原始的fi4_1和fi4_2进行比较

• 模2*pi，th4_1应等于fi4_1等，由np.allclose测试以说明数值舍入误差