סדנאת פייתון לפיזיקאים

דוגמה: מטוטלת כפולה

הפקולטה לפיזיקה, הטכניון. חורף 2013

מרצה: רונן אברבנאל

double penedlum

נרצה לבנות סימולציה של מטוטלת כפולה. לו היינו רוצים לבנות סימולציה של מטוטלת פיזיקאלית פשוטה, זה היה קל, המשוואות, בקירובים סבירים, פתירות אנליטית, והיינו רק צריכים לצייר את העסק.

בגלל שאנחנו מדברים על בעיה שאינה פתירה אנליטית, נצטרך להתאמץ קצת יותר, ולפתור את הבעיה נומרית.

double penedlum

נתחיל בהגדרת הקבועים:

In [3]:

import scipy.constants
g = scipy.constants.g #in SI

החבילה scipy כוללת אוסף חביב של קבועים פיזיקאליים.
לא חייבים, אבל נחמד להשתמש בהם במקום להתחיל להקליד מספרים ארוכים.

נגדיר קבועים נוספים, שרירותיים, שנדרש להם בהמשך:

In [5]:

L1 = 1.0 #in m
L2 = 1.0
M1 = 1.0 #in kg
M2 = 1.0

נגדיר את אינטרוול הזמנים בין הדגימות שלנו לאינטרוול זמן קצר, וניצור מערך למשך 20 שניות, באינטרוול הנדגם,

In [6]:

import numpy as np
dt = 0.05 
t = np.arange(0.0, 20, dt)

נקבע תנאי התחלה: th1, th2 הם הזוויות של שני המוטות, ואילו w1, w2 הן המהירויות הזוותיות.

תנאי ההתחלה הם במעלות/ מעלות לשניה.

In [7]:

th1 = 120.1
w1 = 0.0
th2 = -10.0
w2 = 0.0

נשים בצד משתנה שיעזור לנו להמיר בין מעלות לרדיאנים. גם קבוע כזה יש לנו מוכן.

In [8]:

rad = scipy.constants.degree #=pi/180

נמיר את תנאי ההתחלה לרדינאים, וניצור מערך שכולל את כל תנאי ההתחלה.

In [9]:

state = np.array([th1, w1, th2, w2])*rad
state

Out[9]:

array([ 2.09614043,  0.        , -0.17453293,  0.        ])

לפני שנפתור, נותר לנו החלק הכי חשוב: לתאר את משוואת התנועה.

משוואת דיפרנציאליות רגילות באופן כללי יכולות להכתב בצורה \[]

כאשר במקרה שלנו, \[ ]

עלינו לכתוב פונקציה שמקבלת את $\theta_1,\omega_1,\theta_2,\omega_2,t$ ומחשבת (את הערך המספרי) של ארבעת הפונקציות $f_1,f_2,f_3,f_4$

במקרה של מטוטלת כפולה, הנגזרת של המקיום היא פשוטה למדי, $\frac{{\rm d}\theta_{1}}{{\rm d}t}=f_{1}=\omega_{1}$ הנגזרת עבור $\omega_1$ קצת יותר מסובכת:

$ =f_{2}= $

שתי המשוואות האחרות, דומות למדי.

הקוד:

In [10]:

from numpy import cos, sin, pi
def derivs(state, t):
    
    #unpack values
    (th1, w1, th2, w2) = state
    
    
    dth1_dt = w1

    del_ = th2-th1 
    den1 = (M1+M2)*L1 - M2*L1*cos(del_)*cos(del_)
    dw1_dt = (M2*L1*w1*w1*sin(del_)*cos(del_)
               + M2*g*sin(th2)*cos(del_) + M2*L2*w2*w2*sin(del_)
               - (M1+M2)*g*sin(th1))/den1

    dth2_dt = w2

    den2 = (L2/L1)*den1
    dw2_dt = (-M2*L2*w2*w2*sin(del_)*cos(del_)
               + (M1+M2)*g*sin(th1)*cos(del_)
               - (M1+M2)*L1*w1*w1*sin(del_)
               - (M1+M2)*g*sin(th2))/den2

    #return value 
    dstate_dt = np.array([dth1_dt, dw1_dt, dth2_dt, dw2_dt])
    return dstate_dt

אחרי שהגדרנו את המשוואה הדיפרנציאלית, אפשר לבצע אינטגרציה שלה. נשתמש בפונקציה odeint שמקבלת פונקציית-נגזרות, כמו זו שכתבנו הרגע, מצב התחלתי ומערך המכיל רשימת-זמנים:

In [11]:

import scipy.integrate as integrate
y = integrate.odeint(derivs, state, t)
y

Out[11]:

array([[  2.09614043,   0.        ,  -0.17453293,   0.        ],
       [  2.08363098,  -0.50010351,  -0.1803995 ,  -0.23198823],
       [  2.04618857,  -0.9966524 ,  -0.19719331,  -0.43162774],
       ..., 
       [  0.03002554,   1.52773269,  21.05856252,   5.76938017],
       [  0.12073415,   2.07067146,  21.34396683,   5.64822006],
       [  0.2340227 ,   2.4301546 ,  21.62335263,   5.52660716]])

עכשיו, המערך y הוא מערך דו ממדי: כל אחת מהשוורת מכילה את הפרמטרים $\theta_1,\omega_1,\theta_2,\omega_2,t$ וכל אחת מהעמודות מכילה את הפרמטרים האלו בזמן אחר.

נבצע ציור "נאיווי" של התוצאות:

In [12]:

%pylab inline
import matplotlib.pyplot as plt

plt.plot(t,y)
plt.xlabel("time (secs)")

Populating the interactive namespace from numpy and matplotlib

Out[12]:

<matplotlib.text.Text at 0x44f2e90>

נצייר ציור קצת יותר סביר: הפעם אנחנו רוצים להציג רק את הזוויות:

In [13]:

import matplotlib.pyplot as plt

#We choose only columns 0 and 2
lines = plt.plot(t,y[:,[0,2]])

#labels
plt.xlabel("time (secs)")
plt.ylabel("Angle (radians)")
plt.title("Pendelum angle as function of time")

#Add legend and label lines by order. 
plt.legend(lines, [r"$\theta_1$",r"$\theta_2$"],loc="best")

Out[13]:

<matplotlib.legend.Legend at 0x46e0f10>

תמונה נחמד אחרת תהיה של "מרחב הפאזה" (או לפחות, חלק ממנו)

In [14]:

plt.plot(y[:,0],y[:,2])
plt.xlabel("$\\theta_1$")
plt.ylabel("$\\theta_2$")

Out[14]:

<matplotlib.text.Text at 0x4e1b5d0>

אבל אחרי כל ציורי הביניים הנחמדים האלה, אנחנו רוצים סרטים! אנחנו אוהבים סרטים! סרטים זה כיף!

דבר ראשון, נמיר את הזוויות שיש לנו לערכי x,y של קצוות המוטות של המטוטלת:

In [15]:

x1 = L1 * sin(y[:,0])
y1 = -L1*cos(y[:,0])

x2 = L2 * sin(y[:,2]) + x1
y2 = -L2 * cos(y[:,2]) + y1

double penedlum

נשים לב מה קיבלנו: אם אנחנו רוצים לצייר את המטוטלת באיטרציה i, נצטרך לצייר את:

In [16]:

i = 42
thisx = [0,x1[i],x2[i]]
thisy = [0,y1[i],y2[i]]
plt.plot(thisx,thisy,"-o")
plt.xlim([-2,2])

Out[16]:

(-2, 2)

אבל בואו נעשה אנימציה:

In [17]:

import matplotlib.animation as animation


fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
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)

ani.save('double_pendulum.mp4', fps=15, clear_temp=True)
plt.show()

אבל רגע, זו לא אנימציה! זה סטטי! ה-notebook לא יודע להציג (בנתיים) אנימציה. אבל שמרנו קובץ אז בואו נפתח אותו:

In [18]:

from IPython.display import HTML
video = open("double_pendulum.mp4", "rb").read()
video_encoded = video.encode("base64")
video_tag = '<video controls alt="test" src="data:video/x-m4v;base64,{0}">'.format(video_encoded)
HTML(video_tag)

Out[18]:

האמת היא, שכל העסק הזה הוא נחמד למצגת, אבל הוא הרבה יותר נחמד בתור סקריפט עצמאי:

In [19]:

%load ../example_code/double_pendelum/double.py

In []:

# Double pendulum formula translated from the C code at
# http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c
# Taken from http://matplotlib.org/examples/animation/double_pendulum_animated.html

from numpy import sin, cos, pi, array
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


from numpy import cos, sin, pi
def derivs(state, t):
    
    #unpack values
    (th1, w1, th2, w2) = state
    
    
    dth1_dt = w1
    
    del_ = th2-th1 
    den1 = (M1+M2)*L1 - M2*L1*cos(del_)*cos(del_)
    dw1_dt = (M2*L1*w1*w1*sin(del_)*cos(del_)
               + M2*g*sin(th2)*cos(del_) + M2*L2*w2*w2*sin(del_)
               - (M1+M2)*g*sin(th1))/den1
    
    dth2_dt = w2
    
    den2 = (L2/L1)*den1
    dw2_dt = (-M2*L2*w2*w2*sin(del_)*cos(del_)
               + (M1+M2)*g*sin(th1)*cos(del_)
               - (M1+M2)*L1*w1*w1*sin(del_)
               - (M1+M2)*g*sin(th2))/den2
    
    #return value
    dstate_dt = np.array([dth1_dt, dw1_dt, dth2_dt, dw2_dt])
    return dstate_dt

# create a time array from 0..100 sampled at 0.1 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

rad = pi/180

# initial state
state = np.array([th1, w1, th2, w2])*pi/180.

# 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.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)

#ani.save('double_pendulum.mp4', fps=15, clear_temp=True)
plt.show()

In [20]:

# Double pendulum formula translated from the C code at
# http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c
# Taken from http://matplotlib.org/examples/animation/double_pendulum_animated.html

from numpy import sin, cos, pi, array
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


from numpy import cos, sin, pi
def derivs(state, t):
    
    #unpack values
    (th1, w1, th2, w2) = state
    
    
    dth1_dt = w1
    
    del_ = th2-th1 
    den1 = (M1+M2)*L1 - M2*L1*cos(del_)*cos(del_)
    dw1_dt = (M2*L1*w1*w1*sin(del_)*cos(del_)
               + M2*g*sin(th2)*cos(del_) + M2*L2*w2*w2*sin(del_)
               - (M1+M2)*g*sin(th1))/den1
    
    dth2_dt = w2
    
    den2 = (L2/L1)*den1
    dw2_dt = (-M2*L2*w2*w2*sin(del_)*cos(del_)
               + (M1+M2)*g*sin(th1)*cos(del_)
               - (M1+M2)*L1*w1*w1*sin(del_)
               - (M1+M2)*g*sin(th2))/den2
    
    #return value
    dstate_dt = np.array([dth1_dt, dw1_dt, dth2_dt, dw2_dt])
    return dstate_dt

# create a time array from 0..100 sampled at 0.1 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

rad = pi/180

# initial state
state = np.array([th1, w1, th2, w2])*pi/180.

# 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.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)

#ani.save('double_pendulum.mp4', fps=15, clear_temp=True)
plt.show()

<matplotlib.figure.Figure at 0x4bd12d0>

In []: