assignmet13_Harmonic Motion-short.ipynb

Kernel: SageMath (stable)

import numpy as np
import matplotlib.pyplot as plt
%matplotlib -- inline

בניית הדמיה לחלקיק הנע לאורך קו ישר בהשפעת כוח משתנה- תנועה הרמונית

מטרה:

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

ב. הכרת שיטות לפתרון משוואות דיפרנציאליות מסדר שני

תיאור הבעיה:

גוף מונח על מישור אופקי חלק. בקצה אחד מחובר הגוף לקפיץ שהקבוע שלו k. הקצה השני של הקפיץ מחובר לקיר זקוף. במצב המתואר הקפיץ רפוי, והגוף נמצא בראשית הצירים, כך שהמרחק בין הגוף לקיר הוא האורך הרפוי של הקפיץ. מסיטים את הגוף הצידה לנקודה x0 ומקנים לו מהירות v0. בפעילות תחקרו את העתק הגוף כתלות בזמן ובערכים ההתחלתיים של העתק והמהירות. כמו כן תחקרו ואת השפעת גודל צעד הזמן (dt) על יציבות הפתרון לפי האלגוריתם של אוילר ותכירו אלגוריתמים חדשים.


בניית המודל התיאורטי:

1. בטאו באמצעות k ו-x, את הכוח השקול הפועל על הגוף ברגע שמרחקו מהראשית הוא x, והסבירו מדוע כיוון הכוח הפוך תמיד לכיוון ההעתק (החיכוך זניח).

#F=k*x, and the diraction is always in the opposit diraction the body traveled because the coild is pressured in the diraction the body travled, and then returns the force in the opposit diraction of where it was pressured accourding to Nuton's third law

2. הראו כי משוואת התנועה של הגוף היא : $${d^2x \over dt^2}=-{k\over m}x$$

(הנגזרת השנייה של המקום לפי הזמן פרופורציונית להעתק, כשמקדם הפרופורציה הוא $-{k\over m}$).

#the second divertive of x in t is the acceleratio, which is also equal to F/m (accourding to Nuton's third law), and since F is equal to x*k(accourding to my answer in question 1), the equastion above is right:
#a=F/M
#k*x=F
#a=k*x/m
#d**2x/dt*22=a
#d**2x/dt*22=k*x/m

3. הראו שהפונקציה הבאה פותרת את המשוואה הדיפרנציאלית:$$x(t)=A\cdot cos(\sqrt{ k\over m} t +\delta)$$ כאשר ערכו של הקבוע $\delta$ ניתן על ידי: $$\delta = cos^{-1}({x_0\over A})$$

$\delta$ מכונה הפרש המופע. הראו שקשר זה מתקיים.

x=A*np.cos(np.sqrt(k/m)*t+)#it aint right because you dont know how to shift cos

  File "<ipython-input-4-3cbd6087160d>", line 1
    x=A*np.cos(np.sqrt(k/m)*t+)#it aint right because you dont know how to shift cos
                              ^
SyntaxError: invalid syntax

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

4. פתחו ביטוי למהירות כתלות בזמן.

v=-sqrt(k/m)*A*sin(sqrt(k/m)*t+δ)

5. על מנת למצוא את הפתרון הפרטי יש להתחשב גם במהירות הגוף ההתחלתית. רישמו משוואה נוספת שתקשור בין משתני המצב $A$, $\delta$ לבין $v_0$.

vmax=sqrt(k/m)*A
v0=vmax*sin(sqrt(k/m)*t+δ)=-sqrt(k/m)*A*sin(sqrt(k/m)*t-δ)#probably wrong

  File "<ipython-input-5-0aa9bfd5ac34>", line 2
    v0=vmax*sin(sqrt(k/m)*t+δ)=-sqrt(k/m)*A*sin(sqrt(k/m)*t-δ)#probably wrong
      ^
SyntaxError: can't assign to operator

הערה: הקוד לא עובד

6. נניח כי נתון כי בזמן t=0 הגוף נמצא בנקודה $x_0=0.1 m$ ומהירותו שם אפס. מסת הגוף $m=0.5 kg$ וקבוע הקפיץ $k=15 {N\over m}$. מיצאו את ערכי משתני המצב $A$ ו-$\delta$, ורישמו את הפתרון הפרטי.

A=0.1 because v=0 in that spot, and the amplitude is in the spot where v=0
v=-sqrt(k/m)*A*sin(sqrt(k/m)*t+δ)
A=0.1
v=0
k=15
m=0.5z  
0=-sqrt(15/0.5)*0.1*sin(sqrt(15/0.5)*0+δ)
sin(δ)=0
δ=0

  File "<ipython-input-6-f45ee195c57a>", line 1
    A=0.1 because v=0 in that spot, and the amplitude is in the spot where v=0
                ^
SyntaxError: invalid syntax

הערה: הקוד לא עובד

7. השתמשו ב-numpy על מנת להגדיר את פונקציות המקום והמהירות במקרה הפרטי המתואר בשאלה הקודמת ושרטטו שני גרפים, האחד מתאר את מקום הגוף כתלות בזמן והשני את מהירותו כתלות בזמן. ציירו לפחות שלושה מחזורים.השוו את הגרפים שקיבלתם לגרפים המדוייקים.

import numpy as np
import matplotlib.pyplot as plt
k  =  15#N/m
m =  0.5#kg
x0 = 0.1#m
tf = 5#sec
A=0.1#m
t = np.linspace(0,tf,100)
position = lambda t: A*np.cos(np.sqrt(k/m)*t)
velocity  = lambda t: -np.sqrt(k/m)*A*np.sin(np.sqrt(k/m)*t)
plt.subplot(121)
plt.plot(t,position(t))
plt.xlabel("t[s]")
plt.ylabel("x[m]")
plt.title("x-t graph")
plt.subplot(122)
plt.plot(t,velocity(t))
plt.xlabel("t[s]")
plt.ylabel("v[m]")
plt.title("v-t graph")
#the graphs are pretty simular to the actual graphs, and it is noticeble in them that whenever v=0 x isat his maximum or minimum

Text(0.5,1,'v-t graph')

8. תנועת הגוף היא מחזורית, מה זמן המחזור?

כשניה (על פי המרחק בין הנקודות השוות בגרף האיקס של טי שסירטתי

הערה: נדרשת תשובה מדוייקת

9. תרשים מופע הוא תרשים המתאר את מהירות הגוף כתלות במקום. שרטטו גרף של מהירות הגוף כתלות במקום.
מה צורת הגרף? כיצד תסבירו צורה זו?

import numpy as np
import matplotlib.pyplot as plt
k  =  15#N/m
m =  0.5#kg
x0 = 0.1#m
tf = 5#sec
A=0.1#m
t = np.linspace(0,tf,100)
vlist=[]
plist=[]
position = lambda t: A*np.cos(np.sqrt(k/m)*t)
plist=position(t)
velocity  = lambda t: -np.sqrt(k/m)*A*np.sin(np.sqrt(k/m)*t)
vlist=velocity(t)
fig = plt.figure(figsize = (5,5))
plt.plot(plist,vlist)
plt.xlabel('position',fontsize=12)
plt.ylabel('velocity',fontsize=12)

Text(0,0.5,'velocity')

הערה: חסרה תשובה מילולית

פיתרון נומרי של משוואה דיפרנציאלית מסדר שני

אפשר לראות משוואה דיפרנציאלית מסדר שני כצמד משוואות דיפרנציאליות סימולטניות. למשל אם נגדיר ${dx \over dt}=v$ נוכל להפוך את המשוואה הדיפרנציאלית ${d^2x \over dt^2}=-{k\over m}x$ לצמד משוואות מסדר ראשון: $${dx \over dt }=v$$ $${dv \over dt}=-{k\over m}x$$ את צמד המשוואות האלו ניתן לפתור בשיטות נומריות שונות. בתרגילים שבהמשך תבחנו את הפתרון הנומרי של משוואות התנועה של תנועה הרמונית באמצעות שיטת אוילר.

10. השלימו את הקוד שבהמשך על מנת לפתור בצורה נומרית את משוואות האוסילטור ההרמוני.

import numpy as np
import matplotlib.pyplot as plt

import numpy as np
import matplotlib.pyplot as plt

k  = 15 #N/m
m = 0.5 #kg
x0 = 0.1 #m
v0 = 0

T = 1#[sec] period
omega = 30#rad/sec
tmax = 10#sec

def Euler(tmax,dt):
    N_t = int(tmax/dt)
    t = np.linspace(0, N_t*dt, N_t+1)
    x = np.zeros(N_t+1)
    v = np.zeros(N_t+1)
    # Initial condition
    v[0] = v0
    x[0] = x0
    # Calculate next step
    for n in range(N_t):
        v[n+1] = v[n]+dt*(-omega**2*A*np.cos(omega*t[n]))  #THIS IS WRONG!  Daria
        x[n+1] = x[n]+dt*(-np.sqrt(k/m)*A*np.sin(np.sqrt(k/m)*t[n])) #THIS IS WRONG!  Daria
    return t,x,v

dt=0.1
t,x,v = Euler(tmax,dt)

fig = plt.figure()
l1, l2 = plt.plot(t, x, 'b-', t, x0*np.cos(omega*t), 'r--')
fig.legend((l1, l2), ('numerical', 'exact'), 'upper left')
plt.xlabel('time')
plt.ylabel('position')

Text(0,0.5,'position')

הערה: יש בעיה בקוד. אין כל התאמה בין הגרף המדוייק (האדום) לנומרי (הכחול). הגרף האדום לא נראה כאוסילטור הרמוני.

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

11. שרטטו את האנרגיה הכוללת של המערכת הפונקציה של הזמן. האם הגרף מתנהג כצפוי?

import numpy as np
import matplotlib.pyplot as plt

import numpy as np
import matplotlib.pyplot as plt

A=0.1
k  = 15 #N/m
m = 0.5 #kg
x0 = 0.1 #m
v0 = 0

T = 1#[sec] period
omega = 30#rad/sec
tmax = 10#sec

def Euler(tmax,dt):
    N_t = int(tmax/dt)
    t = np.linspace(0, N_t*dt, N_t+1)
    x = np.zeros(N_t+1)
    v = np.zeros(N_t+1)
    E = np.zeros(N_t+1)
    # Initial condition
    v[0] = v0
    x[0] = x0
    E[0] = (15/2)*0.1**2
    print(E[0])
    # Calculate next step
    for n in range(N_t):
        v[n+1] = v[n]+dt*(-omega**2*A*np.cos(omega*t[n])) #THIS IS WRONG!  Daria
        x[n+1] = x[n]+dt*(-np.sqrt(k/m)*A*np.sin(np.sqrt(k/m)*t[n])) #THIS IS WRONG!  Daria
        E[n+1] = x[n+1]**2*15/2+0.25*v[n+1]**2
    return t,x,v,E
dt=0.1
t,x,v,E = Euler(tmax,dt)
print (E)
print(t)
plt.plot(t,E)
plt.xlabel('time')
plt.ylabel('energy')

0.075
[  7.50000000e-02   2.03250000e+01   4.03459592e-02   1.90641014e+01
   7.82125051e-02   1.65522366e+01   4.87221261e-01   1.31353907e+01
   1.34775956e+00   9.36312829e+00   3.06496377e+00   6.01120027e+00
   5.74358517e+00   3.31206736e+00   9.03926429e+00   1.44942669e+00
   1.27536854e+01   5.34541732e-01   1.62862412e+01   1.28483877e-01
   1.88791373e+01   8.43326764e-03   2.02455628e+01   7.62288071e-02
   2.00260739e+01   4.88779280e-02   1.81681125e+01   1.70654024e-01
   1.52329799e+01   7.72876585e-01   1.16244840e+01   1.93737971e+00
   7.91245256e+00   4.02497100e+00   4.81568727e+00   7.03594973e+00
   2.46774305e+00   1.05305198e+01   9.63687949e-01   1.42337698e+01
   3.24876274e-01   1.74972708e+01   7.11458661e-02   1.95999449e+01
   6.41467627e-03   2.03409966e+01   7.78219967e-02   1.94702314e+01
   7.81414631e-02   1.70732866e+01   3.27220862e-01   1.38068139e+01
   1.17128671e+00   1.01093060e+01   2.67314330e+00   6.54581521e+00
   5.13067067e+00   3.75806067e+00   8.43128980e+00   1.77329183e+00
   1.20449446e+01   6.02812191e-01   1.56359643e+01   1.93210399e-01
   1.85325504e+01   4.49754280e-02   2.00757888e+01   7.69938397e-03
   2.01682705e+01   8.44293170e-02   1.86755915e+01   1.43776853e-01
   1.58141253e+01   5.68285394e-01   1.23156435e+01   1.69898939e+00
   8.62870245e+00   3.55967119e+00   5.29003713e+00   6.36991187e+00
   2.84826120e+00   9.90091706e+00   1.22226686e+00   1.35427746e+01
   3.49297036e-01   1.69187004e+01   1.19305573e-01   1.93592053e+01
   3.54573883e-02   2.02910512e+01   9.23111481e-03   1.97331462e+01
   1.05660191e-01   1.76677927e+01   2.64296494e-01   1.44293783e+01
   9.13735028e-01   1.08013699e+01   2.36852426e+00   7.21753081e+00
   4.59510001e+00]
[  0.    0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1.    1.1
   1.2   1.3   1.4   1.5   1.6   1.7   1.8   1.9   2.    2.1   2.2   2.3
   2.4   2.5   2.6   2.7   2.8   2.9   3.    3.1   3.2   3.3   3.4   3.5
   3.6   3.7   3.8   3.9   4.    4.1   4.2   4.3   4.4   4.5   4.6   4.7
   4.8   4.9   5.    5.1   5.2   5.3   5.4   5.5   5.6   5.7   5.8   5.9
   6.    6.1   6.2   6.3   6.4   6.5   6.6   6.7   6.8   6.9   7.    7.1
   7.2   7.3   7.4   7.5   7.6   7.7   7.8   7.9   8.    8.1   8.2   8.3
   8.4   8.5   8.6   8.7   8.8   8.9   9.    9.1   9.2   9.3   9.4   9.5
   9.6   9.7   9.8   9.9  10. ]

Text(0,0.5,'energy')

הערה: יש בעיה בקוד.ביקשתי את גרף האנרגיה. זה לא נראה נכון..

הערה: חסרה תשובה מילולית.

12. שרטטו תרשים מופע של הפתרון הנומרי. כיצד משתנה התרשים כאשר משנים את dt?

import numpy as np
import matplotlib.pyplot as plt

import numpy as np
import matplotlib.pyplot as plt

A=0.1
k  = 15 #N/m
m = 0.5 #kg
x0 = 0.1 #m
v0 = 0

T = 1#[sec] period
omega = 30#rad/sec
tmax = 10#sec

def Euler(tmax,dt):
    N_t = int(tmax/dt)
    t = np.linspace(0, N_t*dt, N_t+1)
    x = np.zeros(N_t+1)
    v = np.zeros(N_t+1)
    E = np.zeros(N_t+1)
    # Initial condition
    v[0] = v0
    x[0] = x0
    E[0] = (15/2)*0.1**2
    print(E[0])
    # Calculate next step
    for n in range(N_t):
        v[n+1] = v[n]+dt*(-omega**2*A*np.cos(omega*t[n])) #THIS IS WRONG!  Daria
        x[n+1] = x[n]+dt*(-np.sqrt(k/m)*A*np.sin(np.sqrt(k/m)*t[n])) #THIS IS WRONG!  Daria
        E[n+1] = x[n+1]**2*15/2+0.25*v[n+1]**2
    return t,x,v,E
dt=1
while dt<0.000001:
    tmax=dt*100
    T=dt*10
    t,x,v,E = Euler(tmax,dt)
    print (E)
    print(t)
    plt.plot(t,E)
    plt.xlabel('time')
    plt.ylabel('energy')
    dt=dt/10

plt.plot(t,E)

[<matplotlib.lines.Line2D at 0x7f6f928ac400>]

#the lower dt is the more ups and downs the graph has.

הערה: יש בעיה בקוד. עליכם להציג את הגרפים. אני לא רואה אותם. התשובה לא ברורה. מה זה ups ו-downs?. הגרף שאתם מציגים לא טוב, בגלל באג בקוד המקורי, לפני מספר שאלות, ולכן יהיה קשה להסיק לפיו מסקנות הגיוניות.

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

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

הערה: חסרה תשובה מילולית.

Eeuler-Cromer Method

השיטה בה השתמשנו לפתרון המשוואה הדיפרנציאלית מכונה forward Euler method. שיטה זו מבוססת על חישוב הניגזרת עבור הערכים הנוכחיים של משתני המצב, $x_n,v_n$ והערכת משתני ערכיהם של משתני המצב בצעד הבא: $$y_{n+1}=y_n+f(y_{n},t_{n})dt$$

13. השלימו את המשוואות של שיטת אוילר באופן מפורש עבור האוסילטור ההרמוני:

$$x_{n+1}=??$$ $$v_{n+1}=??$$

14. קיימות שיטות נוספות לפתרון משוואות, הנותנות פתרונות מדוייקים יותר. אחת השיטות האלו נקראת שיטת Euler-Cromer. בשיטה זו נעריך קודם את ערכו של אחד ממשתני המצב בצעד הבא, ולאחר מכן נשתמש בערך שמצאנו כדי לחשב את ערכו של המשתנה השני. אפשר לעשות זאת כך, למשל: $$x_{n+1}=x_n + f_1(x_n,v_n,t)$$ $$v_{n+1}=v_n + f_2(x_{n+1},v_n,t)$$ שימו לב שאת הפונקציה השניה מחשבים באמצעות הערך $x_{n+1}$ שחישבנו הרגע. הפונקציה שבהמשך מיישמת שיטה זו.
    השלימו את הפונקציה כך שהיא תיישם את שיטת אוילר-קרומר.

הערה: חסרה תשובה מילולית. בשאלה 13

#13. x(n+1)=x(n)+f(x,t)d**2t
#    v(n+1)=v(n)+f(v,t)d**3t

def Euler_Cromer(tmax,dt):
    N_t = int(tmax/dt)
    t = np.linspace(0, N_t*dt, N_t+1)
    x = np.zeros(N_t+1)
    v = np.zeros(N_t+1)
    # Initial condition
    v[0] = 0
    x[0] = 0.1
    # Calculate next step
    for n in range(N_t):
        x[n+1] = x[n]+dt*(-np.sqrt(k/m)*A*np.sin(np.sqrt(k/m)*t[n])) #THIS IS WRONG!  Daria
        v[n+1] = v[n]+dt*(-omega**2*A*np.cos(omega*t[n])) #THIS IS WRONG!  Daria
    return t,x,v
#this is probably wrong, because i don't know what x and v are referring to (velocity and accelaration maybe?)

הערה: עליך להעזר בניסוח הנוסחאות למעלה, לנסות ולהבין את הקוד, ולהוסיף את הדברים המתאימים..

15. בידקו את התכנסו הפתרון בפונקציה חדשה באמצעות הגרפים של המקום, המהירות ותרשים המופע. מהי מסקנתכם?

#this mathud is more actuarate than the preveious one

הערה: התשובה לא מנומקת, ולמעשה לא מבוססת על דבר..

Second order Euler Method

רשות

בשיטה זו, במקום להעריך את ערכה של הנגזרת רק בנקודה השמאלית של הקטע, נעריך אותה גם בצד הימני של הקטע, ונשתמש בממוצע ביניהן. כדי להעריך את הנגזרת בצד הימני, צריך קודם להעריך את ערכי המשתנים בצעד הבא. לשם כך עושים צעד ניסיון : $$y_{n+1}^* = y_n + f(y_n)dt$$
נשתמש בערך זה כדי להעריך את הנגזרת בצעד הבא, וכך נקבל את הנוסחה של שיטת אוילר מסדר שני: $$y_{n+1} = y_n + \frac{f(y_n)+f(y_{n+1}^*)}{2}dt$$
תרגיל רשות: כיתבו פונקציה חדשה המיישמת את השיטה החדשה עבור האוסילטור ההרמוני.

הרחבה

תנודות של קפיץ דועכות לאחר זמן ארוך מספיק. במודל הפיזיקלי הנפוץ לתיאור תופעה זו מוסיפים איבר מרסן למשוואה. החיכוך פרופורציוני לגודל המהירות ופועל תמיד בניגוד לכיוונה. לכן במשוואה הדיפרנציאלית עבור התאוצה הוא מיוצג על ידי הביטוי $-cv$ כאשר $c$ מייצג את עוצמת הריסון.

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