mirror of
https://github.com/Estom/notes.git
synced 2026-04-14 10:21:08 +08:00
matplotlib & pandas
This commit is contained in:
@@ -0,0 +1,99 @@
|
||||
# 双摆问题
|
||||
|
||||
这个动画说明了双摆问题。
|
||||
|
||||
双摆公式从 http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c 的C代码翻译而来。
|
||||
|
||||
```python
|
||||
from numpy import sin, cos
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.integrate as integrate
|
||||
import matplotlib.animation as animation
|
||||
|
||||
G = 9.8 # acceleration due to gravity, in m/s^2
|
||||
L1 = 1.0 # length of pendulum 1 in m
|
||||
L2 = 1.0 # length of pendulum 2 in m
|
||||
M1 = 1.0 # mass of pendulum 1 in kg
|
||||
M2 = 1.0 # mass of pendulum 2 in kg
|
||||
|
||||
|
||||
def derivs(state, t):
|
||||
|
||||
dydx = np.zeros_like(state)
|
||||
dydx[0] = state[1]
|
||||
|
||||
del_ = state[2] - state[0]
|
||||
den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
|
||||
dydx[1] = (M2*L1*state[1]*state[1]*sin(del_)*cos(del_) +
|
||||
M2*G*sin(state[2])*cos(del_) +
|
||||
M2*L2*state[3]*state[3]*sin(del_) -
|
||||
(M1 + M2)*G*sin(state[0]))/den1
|
||||
|
||||
dydx[2] = state[3]
|
||||
|
||||
den2 = (L2/L1)*den1
|
||||
dydx[3] = (-M2*L2*state[3]*state[3]*sin(del_)*cos(del_) +
|
||||
(M1 + M2)*G*sin(state[0])*cos(del_) -
|
||||
(M1 + M2)*L1*state[1]*state[1]*sin(del_) -
|
||||
(M1 + M2)*G*sin(state[2]))/den2
|
||||
|
||||
return dydx
|
||||
|
||||
# create a time array from 0..100 sampled at 0.05 second steps
|
||||
dt = 0.05
|
||||
t = np.arange(0.0, 20, dt)
|
||||
|
||||
# th1 and th2 are the initial angles (degrees)
|
||||
# w10 and w20 are the initial angular velocities (degrees per second)
|
||||
th1 = 120.0
|
||||
w1 = 0.0
|
||||
th2 = -10.0
|
||||
w2 = 0.0
|
||||
|
||||
# initial state
|
||||
state = np.radians([th1, w1, th2, w2])
|
||||
|
||||
# integrate your ODE using scipy.integrate.
|
||||
y = integrate.odeint(derivs, state, t)
|
||||
|
||||
x1 = L1*sin(y[:, 0])
|
||||
y1 = -L1*cos(y[:, 0])
|
||||
|
||||
x2 = L2*sin(y[:, 2]) + x1
|
||||
y2 = -L2*cos(y[:, 2]) + y1
|
||||
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
|
||||
ax.set_aspect('equal')
|
||||
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_text
|
||||
|
||||
|
||||
def animate(i):
|
||||
thisx = [0, x1[i], x2[i]]
|
||||
thisy = [0, y1[i], y2[i]]
|
||||
|
||||
line.set_data(thisx, thisy)
|
||||
time_text.set_text(time_template % (i*dt))
|
||||
return line, time_text
|
||||
|
||||
ani = animation.FuncAnimation(fig, animate, np.arange(1, len(y)),
|
||||
interval=25, blit=True, init_func=init)
|
||||
|
||||
plt.show()
|
||||
```
|
||||
|
||||
## 下载这个示例
|
||||
|
||||
- [下载python源码: double_pendulum_sgskip.py](https://matplotlib.org/_downloads/double_pendulum_sgskip.py)
|
||||
- [下载Jupyter notebook: double_pendulum_sgskip.ipynb](https://matplotlib.org/_downloads/double_pendulum_sgskip.ipynb)
|
||||
Reference in New Issue
Block a user