assignment9_analytic_soluton-Daria.ipynb

Kernel: SageMath (stable)

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

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

התפרקות רדיוארטיבית

תהליך רדיו אקטיבי הוא תהליך ספונטני שבו גרעין של יסוד מסוים פולט חלקיקי $\alpha$ או $\beta$ והופך לגרעין של יסוד אחר. במקרה שבו הגרעין יפלוט חלקיק $\gamma$ הוא יעבור למצב אנרגטי נמוך יותר מבלי להשתנות. התהליך הוא אקראי, לא ניתן לחזות מראש את הרגע שבו גרעין מסוים יעבור התפרקות רדיואקטיבית. ההסתברות שגרעין יתפרק בפרק זמן קטן $\Delta t$ נמצאת ביחס ישר לפרק זמן זה. פרוש הדבר, שאם ברגע מסוים נתונים $N$ גרעינים וכעבור פרק הזמן $\Delta t$ התפרקו $\Delta N$ גרעינים אזי מתקיים: . $\frac{ΔN}{N}=-\frac{Δt}{τ}$

הגודל $\tau$ הוא גודל האופייני ליסוד המתפרק ומכונה זמן חיים ממוצע. בעזרת קצת אלגברה המשוואה האחרונה תקבל את הצורה: $\frac{ΔN}{Δt}=-\frac{N}{τ}$
המנה ΔN/Δt מתארת את קצב ההתפרקויות הרדיואקטיביות הממוצע בפרק הזמן $\Delta t$. במידה ומשאיפים את $\Delta t$ לאפס מקבלים את קצב ההתפרקות הרגעי: $$\lim_{Δt→0}⁡\frac{ΔN}{Δt}=-\frac{N}{τ}$$ צד שמאל של המשוואה למעלה הוא הנגזרת של N לפי הזמן ולכן: $$\frac{dN}{dt}=-\frac{N}{τ}$$ משוואה אחרונה זו היא משוואה דיפרנציאלית מסדר ראשון. כדי לדעת באופן מפורש מה התלות של N בזמן צריך לפתור משוואה זו. בעבר הכרתם כבר את הפונקציה המעריכית: $y(x)=Ae^{\lambda x}$. המייחד את הפונקציה המעריכית הוא שנגזרתה נמצאת ביחס ישר לפונקציה. גם קצב ההתפרקויות הוא ביחס ישר למספר הגרעינים בזמן. מכאן נניח כי למשוואה הדיפרנציאלית פתרון כללי מהצורה $N(t)=Ae^{\lambda t}$. בזמן t=0 נקבל:
$$N(0)=Ae^{\lambda\cdot 0}=A\cdot 1=A$$ $$A=N_0$$ המקדם A הוא מספר הגרעינים הרדיואקטיביים בהתחלה. כדי לחשב את ערכו של הקבוע λ מציבים את הפונקציה במשוואה הדיפרנציאלית. מקבלים: $$\frac{d}{dt}N_0e^{\lambda t}=-\frac{1}{\tau}N_0e^{\lambda t}$$ $$\lambda N_0e^{\lambda t}=-\frac{1}{\tau}N_0e^{\lambda t}$$ $$\lambda=-\frac{1}{\tau}$$

<br>

הפונקציה sympy.dsolve

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

import sympy as sp
sp.init_printing()
sp.var('t,lambda,tau,N_0')
N=sp.Function('N')
eq=sp.Eq(N(t).diff(t),-N(t)/tau)
eq

$$\frac{d}{d t} N{\left (t \right )} = - \frac{1}{\tau} N{\left (t \right )}$$

sp.dsolve(eq,N(t))

$$N{\left (t \right )} = C_{1} e^{- \frac{t}{\tau}}$$

decay=_.rhs.subs([("C1",N_0)])
decay

$$N_{0} e^{- \frac{t}{\tau}}$$

תרגילים

1. איזה מהגרפים שבתרשים שלמטה מתאר לדעתכם באופן הנכון ביותר את $N(t)$? נמקו!


#A,because its icline is decreesing, and the function is constantly getting smaller

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

2. זמן מחצית החיים של איזוטופ האורניום $U^{232}$ הוא 68.9 שנים.
א. מצאו את זמן החיים הממוצע של $U^{232}$ .
ב. כמה גרם יישארו מ- 1 גרם אחד של $U^{232}$ אחרי פרק זמן השווה לזמן החיים הממוצע?

#2.a.
#e**(-t/T)=0.5
#2.72**(-68.9/T)=o.5
#ln(0.5)=-68.9/T=-0.69
#t=-68.9/-0.69=99.85
#2.b.
#N=n0*e**(-t/T)=n0*e(-1)=1/2.72

3. התרשים שלמטה הוא תרשים כיוונים של המשוואה הדיפרנציאלית $\frac{dN}{dt}=-\frac{N}{τ}$. ידוע כי בזמן 0 מספר החלקיקים הוא N=5. העתיקו את הגרף לתוכנת ציור (למשל paint) ונסו לשרטט על גבי תרשים החיצים את פתרון המשוואה הדיפרנציאלית. העריכו את ערכו של $\tau$.


עוד על תנאי התחלה

פונקציה קדומה לפונקציה $f(x)$ היא פונקציה $F(x)$ המקיימת $\frac{dF}{dx}=f(x)$. הנגזרת של קבוע שווה לאפס, ולכן כל פונקציה מהצורה $G(x)=F(x)+c$ כש- $c$ הוא קבוע, היא פונקציה קדומה של $f(x)$. נסיק מכך שישנה משפחה שלמה של פונקציות הפותרות את המשוואה הדיפרנציאלית (למעשה ראינו כבר באמצעות תרשים הכיוונים שניתן למצוא אין סוף פיתרונות למשוואה הדיפרנציאלית). בדרך כלל לא מעוניינים במשפחה כולה, אלה בפתרון אחד המתאים למקרה פרטי. קביעת פתרון זה מותנית במציאת ערך יחיד ל-$c$. אם יודעים את ערכה של הפונקציה עבור ערך מסוים של המשתנה הבלתי תלוי אפשר למצוא את $c$.
אם עבור $x=x_0$ מתקיים $G(x_0 )=y_0$ נקבל כי: $c=y_0-F(x_0)$.
ככל שבמשוואה הדיפרנציאלית מופיעות נגזרות מסדר גבוה יותר, כך הפתרון יכיל יותר קבועים.

דוגמה: בתנועה שוות תאוצה פותרים משוואה מהצורה $\frac{d^2x}{dt^2} = a$ . משואה זו מציינת שהנגזרת השנייה של המקום קבועה. פתרונה של משוואה זו הוא: $x=x_0+v_0t+\frac{1}{2} at^2$. כפי שאתם רואים, בפונקצית מקום זמן שני קבועים, המקום והמהירות בזמן אפס.

תרגילים
5.במקרה של התפרקות רדיואקטיבית מי הוא קבוע האינטגרציה? ומה משמעותו הפיסיקלית?

#C is the perment, and its physical meaning is the number of particles at the start of the experiment

6.
א. שרטטו תרשים כיוונים עבור המשוואה $\frac{dN}{dt}=-0.2N$.
ב. נניח כי עבור זמן t=0 נתון כי N=9 . העתיקו את תמונת הגרף לתוכנת Paint והוסיפו לו קו המתאר באופן איכותי את הפתרון המתאים. שימרו תמונה זו כקובץ png והעלו קובץ זה לתיקיה שלכם ב-cocalc.

ג. תארו במילים כיצד ישתנה הפתרון אם במקום 0.2 המקדם של N ישתנה ל- 0.3.

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

#6.a.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

dydx=lambda x, y: -0.2*x

x=np.linspace(0.01, 10, 19)
y=np.linspace(0.01, 10, 19)
X,Y=np.meshgrid(x, y)

theta = np.arctan(dydx(X,Y))
vx=np.cos(theta)
vy=np.sin(theta)

plt.figure(figsize=(7, 7))
plt.quiver(X , Y, vx, vy, pivot = 'tail')
plt.title('fig 1')
#6.b.https://imgur.com/a/BJnP2
#6.c.the incline of the graph would be lower
#6.d.t=3.0234(s)

Text(0.5,1,'fig 1')

6. תגובת שרשרת גרעינית

נתונים $N_0$ גרעינים של חומר מסויים (גרעיני אם). החומר הזה דועך בזמן חיים ממוצע $\tau_1$ (קצב דעיכה $\lambda_1=1/\tau_1$) לגרעין יסוד רדיואקטיבי אחר בעל זמן חיים ממוצע $\tau_2$ (קצב דעיכה $\lambda_2=1/\tau_2$). התוצרים של הגרעין השני, גרעין הבת, יציבים. כיצד נמצא את מספר גרעיני הבת כפונקציה של הזמן?

א. הניחו שגרעין הבת איננו רדיו-אקטיבי, כלומר $\lambda_2=0$. כיתבו פונקציה המתארת את מספר גרעיני הבת כפונקציה של הזמן.

#N=N0

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

#dN/dt=0

ג. כעת רישמו את המשוואה הדיפרנציאלית המתארת את קצב השינוי בכמות גרעיני הבת עבור המקרה בו גרעיני הבת גם הם רדיואקטיביים, כלומר $\lambda_2 > 0$.

#dN/dt=1/T2*N

איך נראה הפתרון?

גרעין הרדיום $^{225}_{88}Ra$ עובר התפרקות גרעינית בזמן מחצית חיים של 14.9 ימים והופך לגרעין אקטיניום $^{225}_{89}Ac$.
גרעין האקטיניום הפוך בתורו לגרעין פרנקיום בזמן מחצית חיים של 10 ימים.

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

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

ו. מיצאו פתרון למשוואה הדיפרנציאלית באמצעות dsolve, וציירו אותו על גבי תרשים הכיוונים.

שרשרת ההתפרקויות המלאה שמתחילה בגרעין נפטוניום (Np) מתוארת כאן (https://en.wikipedia.org/wiki/Decay_chain).

#6.d.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
dydx=lambda x, y: -((1/(21.5))*x)

x=np.linspace(0.01, 10, 19)
y=np.linspace(0.01, 10, 19)
X,Y=np.meshgrid(x, y)

theta = np.arctan(dydx(X,Y))
vx=np.cos(theta)
vy=np.sin(theta)

plt.figure(figsize=(7, 7))
plt.quiver(X , Y, vx, vy, pivot = 'tail')
plt.title('fig 1')
#6.e.the amoout of actoniom particles is equal to the derivative of the function of the decay in Radiom particles
#6.f.
T1,N,N0=sp.symbols('T1 N N0')
eq=(N0-(1/T1)*N(T1),N(T1))
answer=sp.dsolve(eq,N(T1))
answer

(N0 - N(T1)/T1, N(T1))

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-2-e1ea43f052e3> in <module>()
     22 eq=(N0-(1/T1)*N(T1),N(T1))
     23 print(eq)
---> 24 answer=sp.dsolve(eq,N(T1))
     25 answer
/usr/local/lib/python3.6/dist-packages/sympy/solvers/ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    583     """
    584     if iterable(eq):
--> 585         match = classify_sysode(eq, func)
    586         eq = match['eq']
    587         order = match['order']
/usr/local/lib/python3.6/dist-packages/sympy/solvers/ode.py in classify_sysode(eq, funcs, **kwargs)
   1492         raise ValueError("classify_sysode() works for systems of ODEs. "
   1493         "For scalar ODEs, classify_ode should be used")
-> 1494     t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
   1495 
   1496     # find all the functions if not given
IndexError: list index out of range

מים פורצים ממכל

נתון מכל גלילי שבאחת מדופנותיו חור. גובהם הרגעי של פני המים במכל יסומן באות D וגובה החריר יסומן ב- h (ראו תרשים). שטח בסיס המכל הוא S ושטח החריר s. אם מניחים שאובדני האנרגיה כתוצאה ממערבולות וחיכוך בחריר זניחים אפשר להראות שמהירות פריצת המים מהחריר, v, מקיימת את המשוואה :$$v=\sqrt{2g(D-h)}$$ (g מציין את תאוצת הנפילה החופשית). משוואה זו מכונה נוסחת טוריצ'לי . בתרגיל שבהמשך תוכיחו נוסחא זו בהסתמך על חוק שימור האנרגיה.


תרגילים

7. במשך פרק זמן קצר $∆t$ פחתה כמות המים במכל בשיעור $\Delta m$. מה הפחת באנרגיה הפוטנציאלית של המים במכל? בטאו את תשובתכם באמצעות $∆m$, $g, h$ ו- $D$.

#7.g*(d-h)*dm

8. נניח שהמהירות בה פורצים המים מהמכל היא $v$, איזה מרחק עברו המים הפורצים מהחריר במשך פרק זמן קצר $∆t$?

#x=v*dt

9. רשמו ביטוי למסת המים $∆m$ שפרצה מהמכל במשך פרק הזמן $∆t$ כתלות בצפיפותם $ρ$ , שטח החריר $s$ , מהירותם $v$ ופרק הזמן הקצר $∆t$.

#dm=dt*v*s

10. בהנחה שכל האנרגיה הפוטנציאלית שאבדה למים במכל הפכה לאנרגיה קינטית של המים הפורצים, הוכיחו את נוסחת טוריצ'לי.

g*(D-h)*dm=Ek
g*h*dm=0.5*dmv**2
v=sqrt(2g(D-h))

11. כאשר מנקבים חריר בדופן האנכית של מכל פלסטיק פורצים המים בכיוון אופקי. באיזה גובה מעל הבסיס צריך לנקב את החריר כדי שלקילוח המים הפורצים יהיה טווח התחלתי מקסימלי (זרם המים יגיע למרחק המקסימלי מדופן המיכל). בטאו את תשובתכם באמצעות גובה פני המים D.

#D-h=gt**2
#h=D-gt**2
#t=sqrt((d-h)/g)
#v=sqrt(2g(D-h))
#v*t=x
#the distance is the biggest when the incline is equal to zero
#d-h=T
#x.diff=1/2*sqrt(T/g)*sqrt(g*T)+sqrt(t/g)*1/2sqrt(g*t)

  File "<ipython-input-35-383f2ac9e8ce>", line 4
    v=sqrt(2g(D-h))
            ^
SyntaxError: invalid syntax

את נפח המים הפורץ מהחריר אפשר לחשב באמצעות הנוסחה: $∆V=sv∆t$. מסמנים את גובה פני המים במכל ב- $D$ ולכן השינוי בגובה פני המים יסומן ב- $∆D$. הפחת בנפח המים שווה לנפח המים שפרצו ולכן $$sv∆t=-S∆D$$לאחר חלוקה ב- $∆t$ והשאפתו לאפס מקבלים את המשואה הדיפרנציאלית: $$\frac{dD}{dt} = -\frac{sv}{S}$$.

אם נציב את נוסחת טוריצ'לי במקום $v$ נקבל: $$\frac{dD}{dt}=-\frac{s}{S}\sqrt{2g(D-h)}$$

גם משוואה זו היא משוואה דיפרנציאלית מסדר ראשון.

תרגילים

12. האם יתכן מצב בו $D$ יהיה קטן מ- $h$?

כן,אבל רק מתי שג׳י שווה לאפס (שאין כוח משיכה בכלל

13. מדוע מופיע סימן המינוס במשוואה הדיפרנציאלית?

#because the surface is getting lower, since water only come out and not in

14. נתון מכל מים ששטח בסיסו $S=400cm^2$ ובתחתיתו חריר ששטחו $s=1cm^2$. בזמן $t=0$ גובה המים במכל $D_0=1m$. מצאו ביטוי לגובה המים במכל כתלות בזמן. שרטטו גרף של הגובה כתלות בזמן וחשבו כעבור כמה זמן המכל יתרוקן.

import sympy as sp
import matplotlib.pyplot as pls
sp.var('D,lambda,t,s,S,D_0')
D=sp.Function('D')
D_0=1
eq=sp.Eq(D(t).diff(t),(s/S)*sp.sqrt(2*g*D)
t_list=range(1,100)
D_list=[D_0]
g=10
for i in t_list:
    eq.subs([("D",D_list[-1]),(t,i)])
    new_D=sp.dsolve(eq,D(t))
    D_list.append(new_D)
pls.plot(D_list,t_list)

  File "<ipython-input-24-8c8c230c24ca>", line 7
    t_list=range(1,100)
         ^
SyntaxError: invalid syntax

דוגמאות נוספות

16. משך הזמן הדרוש לביצוע פעולה מסוימת הולך ומתקצר ככל שמתאמנים בביצוע המשימה. ככל שמספר האימונים גדל, זמן הביצוע קטן עד ששהוא מתייצב על ערך מסוים $M$. נסמן את מספר האימונים ב- $N$ ואת זמן הביצוע כ- $P(N)$ (זמן הביצוע הוא פונקציה של מספר האימונים $N$). המשוואה הדיפרנציאלית $\frac{dP(N)}{dN}=k·(M-P(N))$ מתארת באופן מקורב את תלות קצב השתנות זמן הביצוע במספר האימונים. $k$ הוא קבוע מסוים התלוי בסוג המיומנות הנלמדת. הסבירו את המשוואה הדיפרנציאלית. בעזרת sympy מצאו את פתרון המשוואה. הזמן הממוצע הדרוש לביצוע משימה מסוימת לאחר ניסיון ראשון הוא 4 דקות, לאחר הניסיון השני 3.7 דקות ולאחר ניסיון שלישי 3.5 דקות. שרטטו עקומת למידה לביצוע משימה זו. מצאו מה הזמן הקטן ביותר אליו ניתן לשאוף בביצוע משימה זו.

17. מודל Gompertz הוא אחד המודלים המוצלחים לחיזוי של קצב הגידול של גידול סרטני. המודל מבוסס על המשוואה הדיפרנציאלית: $\frac{dN}{dt}=k\cdot N\cdot ln⁡(\frac{M}{N})$ כאשר $N(0)=N_0$. מספר תאי הסרטן מסומן ב-$N$, הפרמטרים $k$ ו-$M$ הם קבועים.
א. אם הזמן נימדד בדקות באיזה יחידות נמדד $k$?
ב. ככל ש- $N$ שואף ל-M מה קורה לקצב הגידול?
ג. מצאו בעזרת sympy פיתרון למשוואת Gompertz. ד. שרטטו גרף של מספר תאי הסרטן N כתלות בזמן עבור מצב שבו $k=1$ , $N(0)=30$ ו- $M=109$.