微积分基础¶
微积分是理解机器人运动和控制的数学基础。本章将讲解机器人学所需的微积分知识。
学习目标¶
完成本章后，你将能够：
- 理解导数和微分的概念
- 掌握积分的基本方法
- 理解微分方程的求解
- 应用微积分解决机器人学问题
1. 导数基础¶
1.1 导数的定义¶
导数描述了函数在某一点的变化率：
1.2 基本导数公式¶
import sympy as sp
x = sp.Symbol('x')
# 基本导数
f1 = x**2
f1_diff = sp.diff(f1, x) # 2*x
f2 = sp.sin(x)
f2_diff = sp.diff(f2, x) # cos(x)
f3 = sp.exp(x)
f3_diff = sp.diff(f3, x) # exp(x)
f4 = sp.log(x)
f4_diff = sp.diff(f4, x) # 1/x
print("x² 的导数:", f1_diff)
print("sin(x) 的导数:", f2_diff)
print("exp(x) 的导数:", f3_diff)
print("log(x) 的导数:", f4_diff)
1.3 导数规则¶
import sympy as sp
x = sp.Symbol('x')
f = sp.Function('f')(x)
g = sp.Function('g')(x)
# 和法则
# (f + g)' = f' + g'
sum_rule = sp.diff(f + g, x)
# 积法则
# (f * g)' = f' * g + f * g'
product_rule = sp.diff(f * g, x)
# 商法则
# (f / g)' = (f' * g - f * g') / g²
quotient_rule = sp.diff(f / g, x)
# 链式法则
# [f(g(x))]' = f'(g(x)) * g'(x)
chain_rule = sp.diff(f.subs(x, g), x)
print("和法则:", sum_rule)
print("积法则:", product_rule)
print("商法则:", quotient_rule)
1.4 偏导数¶
import sympy as sp
x, y = sp.symbols('x y')
# 多元函数
f = x**2 * y + sp.sin(x * y)
# 对 x 的偏导数
df_dx = sp.diff(f, x) # 2*x*y + y*cos(x*y)
# 对 y 的偏导数
df_dy = sp.diff(f, y) # x**2 + x*cos(x*y)
# 二阶偏导数
d2f_dx2 = sp.diff(f, x, 2) # 2*y - y**2*sin(x*y)
d2f_dxdy = sp.diff(f, x, y) # 2*x + cos(x*y) - x*y*sin(x*y)
print("∂f/∂x:", df_dx)
print("∂f/∂y:", df_dy)
print("∂²f/∂x²:", d2f_dx2)
print("∂²f/∂x∂y:", d2f_dxdy)
2. 梯度与方向导数¶
2.1 梯度¶
梯度是偏导数组成的向量，指向函数增长最快的方向：
import numpy as np
def gradient(f, x, h=1e-5):
"""数值梯度"""
n = len(x)
grad = np.zeros(n)
for i in range(n):
x_plus = x.copy()
x_minus = x.copy()
x_plus[i] += h
x_minus[i] -= h
grad[i] = (f(x_plus) - f(x_minus)) / (2 * h)
return grad
# 示例：f(x, y) = x² + y²
def f(x):
return x[0]**2 + x[1]**2
# 在点 (1, 2) 处的梯度
x = np.array([1.0, 2.0])
grad = gradient(f, x)
print("梯度:", grad) # [2, 4]
# 梯度方向是函数增长最快的方向
# 负梯度方向是函数下降最快的方向
2.2 方向导数¶
方向导数描述了函数在某个方向上的变化率：
import numpy as np
def directional_derivative(f, x, direction):
"""计算方向导数"""
grad = gradient(f, x)
# 单位化方向向量
u = direction / np.linalg.norm(direction)
return np.dot(grad, u)
# 示例
def f(x):
return x[0]**2 + x[1]**2
x = np.array([1.0, 2.0])
# 不同方向的方向导数
directions = [
np.array([1, 0]), # x 方向
np.array([0, 1]), # y 方向
np.array([1, 1]), # 45度方向
]
for d in directions:
dd = directional_derivative(f, x, d)
print(f"方向 {d}: 方向导数 = {dd:.2f}")
3. 积分基础¶
3.1 不定积分¶
不定积分是导数的逆运算：
import sympy as sp
x = sp.Symbol('x')
# 基本积分
f1 = x**2
F1 = sp.integrate(f1, x) # x**3/3
f2 = sp.sin(x)
F2 = sp.integrate(f2, x) # -cos(x)
f3 = sp.exp(x)
F3 = sp.integrate(f3, x) # exp(x)
f4 = 1/x
F4 = sp.integrate(f4, x) # log(x)
print("∫x²dx =", F1)
print("∫sin(x)dx =", F2)
print("∫exp(x)dx =", F3)
print("∫1/x dx =", F4)
3.2 定积分¶
定积分计算函数在区间上的累积量：
import sympy as sp
x = sp.Symbol('x')
# 定积分
f = x**2
area = sp.integrate(f, (x, 0, 1)) # 1/3
print("∫[0,1] x²dx =", area)
# 数值积分
import numpy as np
from scipy import integrate
def f(x):
return x**2
result, error = integrate.quad(f, 0, 1)
print("数值积分结果:", result)
print("误差估计:", error)
3.3 多重积分¶
import sympy as sp
x, y = sp.symbols('x y')
# 二重积分
f = x * y
area = sp.integrate(f, (x, 0, 1), (y, 0, 1)) # 1/4
print("∫∫ xy dxdy =", area)
# 数值二重积分
import numpy as np
from scipy import integrate
def f(y, x):
return x * y
result, error = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1)
print("数值二重积分:", result)
4. 微分方程¶
4.1 常微分方程 (ODE)¶
import sympy as sp
# 符号求解
x = sp.Symbol('x')
y = sp.Function('y')
# 一阶 ODE: y' = y
eq1 = sp.Eq(y(x).diff(x), y(x))
sol1 = sp.dsolve(eq1, y(x))
print("y' = y 的解:", sol1)
# 二阶 ODE: y'' + y = 0
eq2 = sp.Eq(y(x).diff(x, 2) + y(x), 0)
sol2 = sp.dsolve(eq2, y(x))
print("y'' + y = 0 的解:", sol2)
4.2 数值求解 ODE¶
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# 定义 ODE: y' = -2y
def model(y, t):
dydt = -2 * y
return dydt
# 初始条件
y0 = 1
# 时间点
t = np.linspace(0, 5, 100)
# 求解 ODE
y = odeint(model, y0, t)
# 绘制结果
plt.plot(t, y)
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.title("Solution of y' = -2y")
plt.grid(True)
plt.show()
4.3 机器人动力学中的微分方程¶
import numpy as np
from scipy.integrate import odeint
# 单摆动力学
# θ'' + (g/l)sin(θ) = 0
def pendulum(state, t, g=9.81, l=1.0):
theta, omega = state
dtheta_dt = omega
domega_dt = -(g/l) * np.sin(theta)
return [dtheta_dt, domega_dt]
# 初始条件：角度=30度，角速度=0
theta0 = np.radians(30)
omega0 = 0
state0 = [theta0, omega0]
# 时间点
t = np.linspace(0, 10, 1000)
# 求解
solution = odeint(pendulum, state0, t)
# 绘制结果
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(t, np.degrees(solution[:, 0]))
plt.xlabel('Time (s)')
plt.ylabel('Angle (degrees)')
plt.title('Pendulum Angle')
plt.grid(True)
plt.subplot(1, 2, 2)
plt.plot(t, solution[:, 1])
plt.xlabel('Time (s)')
plt.ylabel('Angular velocity (rad/s)')
plt.title('Pendulum Angular Velocity')
plt.grid(True)
plt.tight_layout()
plt.show()
5. 泰勒展开¶
5.1 泰勒级数¶
import sympy as sp
x = sp.Symbol('x')
a = 0 # 展开点
# 泰勒展开
f = sp.sin(x)
taylor = sp.series(f, x, a, n=6) # 展开到 6 阶
print("sin(x) 的泰勒展开:", taylor)
# 数值验证
import numpy as np
def taylor_sin(x, n=6):
"""sin(x) 的泰勒近似"""
result = 0
for i in range(n):
term = ((-1)**i) * (x**(2*i+1)) / np.math.factorial(2*i+1)
result += term
return result
# 比较
x_val = 0.5
print(f"sin({x_val}) = {np.sin(x_val):.6f}")
print(f"泰勒近似 = {taylor_sin(x_val):.6f}")
print(f"误差 = {abs(np.sin(x_val) - taylor_sin(x_val)):.2e}")
5.2 应用：线性化¶
import numpy as np
def linearize(f, x0, h=1e-5):
"""函数在 x0 处的线性化"""
f0 = f(x0)
df = (f(x0 + h) - f(x0 - h)) / (2 * h)
def linear_approx(x):
return f0 + df * (x - x0)
return linear_approx
# 示例：线性化 sin(x) 在 x=0 处
f = np.sin
x0 = 0
linear_f = linearize(f, x0)
# 比较
x = np.linspace(-1, 1, 100)
plt.plot(x, f(x), label='sin(x)')
plt.plot(x, linear_f(x), label='Linear approximation')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.title('Linearization of sin(x) at x=0')
plt.show()
6. 应用示例¶
6.1 机器人速度分析¶
import numpy as np
def forward_kinematics(theta, l1=1.0, l2=0.8):
"""正运动学"""
theta1, theta2 = theta
x = l1 * np.cos(theta1) + l2 * np.cos(theta1 + theta2)
y = l1 * np.sin(theta1) + l2 * np.sin(theta1 + theta2)
return np.array([x, y])
def jacobian(theta, l1=1.0, l2=0.8):
"""雅可比矩阵"""
theta1, theta2 = theta
J = np.array([
[-l1*np.sin(theta1) - l2*np.sin(theta1+theta2),
-l2*np.sin(theta1+theta2)],
[l1*np.cos(theta1) + l2*np.cos(theta1+theta2),
l2*np.cos(theta1+theta2)]
])
return J
# 关节角度
theta = np.array([np.pi/4, np.pi/6])
# 关节速度
dtheta = np.array([0.1, 0.2])
# 计算末端速度
J = jacobian(theta)
dx = J @ dtheta
print("关节速度:", dtheta)
print("末端速度:", dx)
6.2 轨迹规划¶
import numpy as np
import matplotlib.pyplot as plt
def cubic_trajectory(t0, tf, q0, qf, v0=0, vf=0):
"""三次多项式轨迹规划"""
# 系数矩阵
A = np.array([
[1, t0, t0**2, t0**3],
[0, 1, 2*t0, 3*t0**2],
[1, tf, tf**2, tf**3],
[0, 1, 2*tf, 3*tf**2]
])
# 目标向量
b = np.array([q0, v0, qf, vf])
# 求解系数
coeffs = np.linalg.solve(A, b)
return coeffs
def evaluate_trajectory(coeffs, t):
"""计算轨迹值"""
a0, a1, a2, a3 = coeffs
q = a0 + a1*t + a2*t**2 + a3*t**3
dq = a1 + 2*a2*t + 3*a3*t**2
ddq = 2*a2 + 6*a3*t
return q, dq, ddq
# 轨迹参数
t0, tf = 0, 2 # 时间
q0, qf = 0, 1 # 位置
v0, vf = 0, 0 # 速度
# 计算系数
coeffs = cubic_trajectory(t0, tf, q0, qf, v0, vf)
# 评估轨迹
t = np.linspace(t0, tf, 100)
q, dq, ddq = evaluate_trajectory(coeffs, t)
# 绘制结果
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.plot(t, q)
plt.xlabel('Time (s)')
plt.ylabel('Position')
plt.title('Position')
plt.grid(True)
plt.subplot(1, 3, 2)
plt.plot(t, dq)
plt.xlabel('Time (s)')
plt.ylabel('Velocity')
plt.title('Velocity')
plt.grid(True)
plt.subplot(1, 3, 3)
plt.plot(t, ddq)
plt.xlabel('Time (s)')
plt.ylabel('Acceleration')
plt.title('Acceleration')
plt.grid(True)
plt.tight_layout()
plt.show()
7. 练习题¶
基础练习¶
- 计算函数 f(x) = x³ + 2x² - 5x + 1 的导数
- 计算定积分 ∫[0,π] sin(x)dx
- 求解微分方程 y' = 2y, y(0) = 1
进阶练习¶
- 计算二元函数 f(x,y) = x²y + xy² 的梯度
- 使用泰勒展开近似 cos(x) 在 x=0 附近
- 实现单摆动力学的数值仿真
应用练习¶
- 实现 2 连杆机器人的速度分析
- 使用三次多项式规划关节轨迹
- 实现简单的轨迹跟踪控制器
