动图之一个圆绕另一个圆转动

已知圆1与圆2，圆2不动，圆1绕圆2转一圈，圆1自转多少圈？

答：圆1的圆心走了 2 π ( r 1 + r 2 ) 2\pi(r_1+r_2) 2π(r1+r2)，圆1自转了的路程也是 2 π ( r 1 + r 2 ) 2\pi(r_1+r_2) 2π(r1+r2)，因此圆1自转了 2 π ( r 1 + r 2 ) 2 π r 1 \frac{2\pi(r_1+r_2)}{2\pi r_1} 2πr12π(r1+r2)圈

动图

# -*- coding: utf-8 -*-
"""
Created on Mon Jun 14 20:22:04 2021@author: Leslie Lee
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animationdef coord(time, step, W, w, r1, r2):'''r1为动圆半径 r2为静圆半径 W为公转角速度 w为自转角速度 time为离散时间长度'''t = np.arange(0, time, step)x = np.sin(W*t)*r1y = np.cos(W*t)*r1x1 = np.sin(W*t)*(r1+r2)y1 = np.cos(W*t)*(r1+r2)x2 = x1 + r1*np.sin(w*t/2)y2 = y1 - r1*np.cos(w*t/2)return x, y, x1, y1, x2, y2time = 500
step = 0.05
r1 = 2
r2 = 2
W = np.deg2rad(3.6)*10
w =np.deg2rad(3.6)*10
x,y,x1,y1,x2,y2 = coord(time, step, W, w, r1, r2)# 建立fig
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-8, 8), ylim=(-8, 8))
ax.grid()
line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs' # 格式化输出 保留一位浮点数
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)def init():line.set_data([], [])time_text.set_text('')return line, time_textdef animate(i):thisx = [0, x[i], x1[i], x2[i]]thisy = [0, y[i], y1[i], y2[i]]line.set_data(thisx, thisy)time_text.set_text(time_template % (i*step))return line, time_textani = animation.FuncAnimation(fig, animate, np.arange(0, len(y)),interval=25, blit=True, init_func=init)
plt.show()

在动画的基础上又将轨迹画了上去

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animationdef coord(time, step, W, w, r1, r2):'''r1为动圆半径 r2为静圆半径 W为公转角速度 w为自转角速度 time为离散时间长度'''t = np.arange(0, time, step)x = np.sin(W*t)*r1y = np.cos(W*t)*r1x1 = np.sin(W*t)*(r1+r2)y1 = np.cos(W*t)*(r1+r2)x2 = x1 + r1*np.sin(w*t/2)y2 = y1 - r1*np.cos(w*t/2)return x, y, x1, y1, x2, y2time = 500
step = 0.05
r1 = 2
r2 = 2
W = np.deg2rad(3.6)*10
w =np.deg2rad(3.6)*10
x,y,x1,y1,x2,y2 = coord(time, step, W, w, r1, r2)fig = plt.figure()
ax = plt.gca()
ax.set_aspect(1)plt.plot(x, y)
plt.plot(x1, y1)
plt.plot(x2, y2)line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs' # 格式化输出 保留一位浮点数
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)def init():line.set_data([], [])time_text.set_text('')return line, time_textdef animate(i):thisx = [0, x[i], x1[i], x2[i]]thisy = [0, y[i], y1[i], y2[i]]line.set_data(thisx, thisy)time_text.set_text(time_template % (i*step))return line, time_textani = animation.FuncAnimation(fig, animate, np.arange(0, len(y)),interval=25, blit=True, init_func=init)
plt.show()

保存gif，需要安装pillow

ani.save('F:\double_pendulum.gif', writer='pillow')

验证

如何验证圆1自转了多少圈？
即验证蓝线上的红点走够一圈，绿线的长度是否为 2 π ( r 1 + r 2 ) 2\pi(r_1+r_2) 2π(r1+r2)
方1. 曲线积分，已知曲线参数方程 ( x 2 ( t ) , y 2 ( t ) ) (x_2(t),y_2(t)) (x2(t),y2(t))，来求弧长
方2. 已知绿线的离散点，对相邻两个离散点求距离，将所有距离累计起来便是弧长
方1出现了不可积项，因此放弃

方1积分程序

# -*- coding: utf-8 -*-
"""
Created on Tue Jun 15 08:44:13 2021@author: Leslie Lee  用python的库 sympy 求积分 https://blog.csdn.net/t4ngw/article/details/105770161
https://wenku.baidu.com/view/5b46565491c69ec3d5bbfd0a79563c1ec4dad71c.html
"""
import numpy as np
from sympy import var, diff, integrate, sin, cos, sqrt, simplify, trigsimp, pidef arc(r1, r2, W, w):var("t")x1 = sin(W*t)*(r1+r2)y1 = cos(W*t)*(r1+r2)x2 = x1 + r1*sin(w*t/2)y2 = y1 - r1*cos(w*t/2)# 求导dx2 = diff(x2, t)dy2 = diff(y2, t)# 积分end = 2*pi/wres = integrate(sqrt(dx2**2+dy2**2), (t, 0, end))res = trigsimp(simplify(res))return resW = np.deg2rad(3.6)*10
w =np.deg2rad(3.6)*10
r1 = 2
r2 = 2
res = arc(2, 2, W, w)

方2累加程序

# -*- coding: utf-8 -*-
"""
Created on Mon Jun 21 11:27:28 2021@author: Leslie Lee
"""import numpy as np
import matplotlib.pyplot as pltdef arc(r1, r2, W, w, num):time = 2*np.pi/Wstep = time/num# print(time, step)t = np.arange(0, time, step)x1 = np.sin(W*t)*(r1+r2)y1 = np.cos(W*t)*(r1+r2)x2 = x1 + r1*np.sin(w*t/2)y2 = y1 - r1*np.cos(w*t/2)# 实际值res = 0for i in range(int(len(x2)-1)):res = res + ((x2[i+1] - x2[i])**2 + (y2[i+1] - y2[i])**2)**0.5# 目标值target_res = 2*np.pi*(r1+r2)return res, target_res# 结果与转速无关 与离散时间的间隔有关
W = np.deg2rad(3.6)*10
w =np.deg2rad(3.6)*10
r1 = 2
r2 = 2nums = []
results = []
target_results = []
for i in range(1, 1000, 10): res, t_res = arc(2, 2, W, w, i)nums.append(i)results.append(res)target_results.append(t_res)
plt.plot(nums, results, label='real')
plt.plot(nums, target_results, label='target')
plt.xlabel('the length of time series')
plt.legend()
plt.show()

验证了确实对了，但时间序列间隔也不是越小越好，有个限度，比如这个例子中间隔为50结果就是正常结果(两条曲线交点)

参考：
求解双摆与单摆运动轨迹并绘制动图。https://blog.csdn.net/qq_37083038/article/details/117884311
python学习之matplotlib绘制动图（FuncAnimation()参数）。https://www.cnblogs.com/zhouzhe-blog/p/9614360.html