Contact Us!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download
Views: 209
Image: ubuntu2004
Tweet
Share
Share
Kernel: Python 3 (Anaconda 2020)
import matplotlib.pyplot as plt
import numpy as np

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

Shlomo Rozenfeld

 אוילר מסדר ראשון

לכל המשוואות שהובאו כאן ישנו פתרון אנליטי. בפתרון אנליטי הכוונה היא שניתן לרשום במפורש כיצד y תלוי ב- x. למשל בהתפרקות רדיואקטיבית אך לא תמיד המצב כך. ישנן משואות דיפרנציאליות שאין אנחנו יודעים לבטא את פתרונן באופן אלגברי. במקרים אלו הפתרון נעשה בדרך נומרית. מוצאים ערך מקורב למשתנה התלוי עבור סידרת ערכים בדידים של המשתנה הבלתי תלוי. נניח שנתונה המשוואה הדיפרנציאלית . ידוע גם ערכו ההתחלתי של המשתנה התלוי y0y_0 עבור ערך מסוים x0x_0, של המשתנה התלוי. רוצים לחשב את ערכו של המשתנה התלוי כאשר ערכו של המשתנה הבלתי תלוי הוא x. השיטה הפשוטה ביותר לעשות זאת היא קרוב אוילר מסדר ראשון . כאשר ידוע כי ערכו של המשתנה התלוי הוא y0y_0 כשערכו של המשתנה הבלתי הוא x0x_0, אם נדע את שיפוע המייתר העובר בין הנקודות (x0,y0)(x_0,y0) ו- (x,y)(x,y) נוכל לחשב את ערך המשתנה התלוי y עבור ערך ידוע של x. הקו המרוסק בתרשים מתאר את הפונקציה האמתית והמיתר המחבר את הנקודה (x0,y0)(x_0,y_0) ונקודה הקרובה אליה הנקודה . ידוע כי שיפוע המייתר הוא s.

pic_7

מהגדרת שיפוע המיתר [removed]s=y−y0x−x0s=\frac{y-y_0}{x-x_0} נקבל: [removed]y=y0+s⋅(x−x0)y=y_0+s\cdot (x-x_0). הבעיה שאין אנו יודעים את שיפוע המיתר. אבל ידועה לנו המשוואה הדיפרנציאלית: [removed]dydx=f(x,y)\frac{dy}{dx}=f(x,y). אם נציב ב- f את x0x_0 ואת y0y_0 נקבל את שיפוע המשיק לגרף בנקודה זו. ככל ש- x יהיה קרוב יותר ל- x0x_0 כך שיפוע המשיק יתקרב לשיפוע המייתר המחבר את הנקודות (x0,y0)(x_0,y_0) ו- (x,y). נוכל לפיכך לכתוב: [removed]y=y0+f(x0,y0)(x−x0)y=y_0+f(x_0,y_0)(x-x_0) נסמן את [removed]x−x0x-x_0 ב-[removed]Δx\Delta x נקבל:[removed]y=y0+f(x0,y0)Δx y=y_0+f(x_0,y_0)\Delta x מהערך החדש של y אפשר לחשב את ערך הפונקציה עבור נקודה במרחק Δx\Delta x וכן הלאה.

pic_8

התרשים שלמעלה ממחיש את התהליך. בהתחלה קובעים את גודל הצעד של המשתנה הבלתי תלוי: [removed]h=Δxh=\Delta x ואת שיעורי נקודת ההתחלה a ו-b בכל פעם מחשבים את ערכו החדש של x ואעת ערכו החדש של y. הנוסחא הכללית של קרוב אוילר היא הנוסחא הבאה: [removed]yn+1=yn+f(xn,yn)⋅hy_{n+1}=y_n+f(x_n,y_n)\cdot h בכל פעם נעשה שימוש לחישוב הערכים החדשים, בערכים שהתקבלו קודם. התהליך מודגם בתרשים:

pic_9

קטע הקוד שבהמשך מתארת את החישוב עבור משוואה דיפרנציאלית מהצורה [removed]dydx=x2\frac{dy}{dx}=x^2. ידוע כי עבור x=3 ערכו של y הוא 1.0 . הבחירה של [removed]Δx=h\Delta x=h תלויה בדיוק הרצוי. ככל שנבחר Δx\Delta x קטן יותר, כך החישוב יהיה מדויק יותר.



#ערכי ההתחלה
x0=3.0
y0=1.0
dx=0.1 #הגדרת גודל הצעד
dydx=lambda x:x*x #חישוב הנגזרת בנקודה מסוימת

xn=x0
yn=y0
print( "| n |  xn   |   yn  |")
print ("|---|-------|-------|")
print( "| 0 | {0:1.2f}  |  {1:1.2f} |".format(xn,yn))
for i in range(1,5):# 
    yn=yn+dydx(xn)*dx
    xn=xn+dx
    print ("| {0} | {1:1.2f}  |  {2:1.2f} |".format(i,xn,yn))
| n | xn | yn | |---|-------|-------| | 0 | 3.00 | 1.00 | | 1 | 3.10 | 1.90 | | 2 | 3.20 | 2.86 | | 3 | 3.30 | 3.89 | | 4 | 3.40 | 4.97 |


1. למשוואה הדיפרנציאלית [removed]dydx=x2\frac{dy}{dx}=x^2 עבור תנאי ההתחלה y(3)=1 פיתרון אנאליטי: [removed]y(x)=x33−8y(x)=\frac{x^3}{3}-8. שנו את הקוד בדוגמה כך שיוסיף לטבלה עמודה של הערך המדויק של y ועמודה של הערך המוחלט של ההפרש בין הערך המדויק לזה המתקבל משיטת אוילר. בדקו איך השגיאה משתנה כאשר משנים את Δx. .


2.הפונקציה שבהמשך פותרת בשיטת אוילר כל משוואה דיפרנציאלית מהצורה. :

א. רשמו הסבר במקומות המסומנים

ב. בדקו את הקוד על הדוגמא שלמעלה

def Euler(f,x0,y0,x1,dx = 0.1):
    x = x0#?
    y = y0#?
    n = int((x1 - x0) / dx +1)#?
    x_array = np.zeros(n)#?
    x_array[0] = x0
    y_array = np.zeros(n)
    y_array[0] = y0
    
    for i in range(1, n):
        y_array[i] = (y_array[i-1]+f(x_array[i-1],y_array[i-1])*dx)#?
        x_array[i] =  x_array[i-1] + dx
    return x_array, y_array


3. נתונה המשוואה הדיפרנציאלית [removed]dNdt=−5N\frac{dN}{dt}=-5N. כמו כן ידוע כי עבור זמן t=0, N(0)=6.0.

א. שרטטו גרף של פתרון נומרי של המשוואה עבור? [removed]0≤t≤0.50\le t\le 0.5 הניחו כי dt=0.1.

ב. מצאו את הפתרון האנליטי למשוואה והוסיפו אותו על הגרף הקודם
.
 ג. צרו גרף חדש של הפתרון הנומרי אבל הפעם עבור dt=0.01. גם הפעם הוסיפו את הפתרון האנליטי לגרף. מה המסקנה?


4. א. בעזרת התרגיל הקודם מצאו את השגיאה היחסית בחישוב N בזמן t=0.5 עבור כל אחד אחד מערכי dt הבאים: 0.01,0.02,0.025,0.05,0.1 .
ב. שרטטו גרף של השגיאה היחסית כתלות ב- dt. מה המסקנה?
הערה: השגיאה היחסית מוגדרת כ- ∣Nnumeric−Ntrue∣Ntrue\frac{|N_{numeric}-N_{true}|}{N_{true}}

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

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

pic


את נפח המים הפורץ מהחריר אפשר לחשב באמצעות הנוסחה: . גובה פני המים במכל סימנו DD ולכן השינוי בגובה פני המים יסומן ב- ∆D∆D. הפחת בנפח המים שווה לנפח המים שפרצו ולכן . לאחר חלוקה ב- ∆t∆t והשאפתו לאפס מקבלים את המשואה הדיפרנציאלית: .

שימוש בנוסחת טירוצ'לי עבור vv נותן: dDdt=−sS2g(D−h)\frac{dD}{dt}=-\frac{s}{S}\sqrt{2g(D-h)}

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

תרגילים

 5. האם יתכן מצב ש- DD יהיה קטן מ- hh?


6 .מדוע מופיע סימן המינוס במשוואות השונות?


7. נתון מכל שקוטר בסיסו 15cm, קוטר החריר 0.5mm והוא בגובה של 12cm מעל לבסיס המכל. בעזרת שיטת אוילר שרטטו גרף את מהירות פריצת המים כתלות בזמן עד הרגע שבו גובה פני המים הוא 0.5m מעל הבסיס. מצאו מה צריך להיות גודל הפסיעה כדי לקבל שגודל השגיאה היחסית לא יעלה על 0.1%.


8. גוף נופל ממנוחה. מסת הגוף 0.5kg קבוע התנגדות האוויר c=0.1kg/m ותאוצת הנפילה החופשית g=9.8m/s^2.הרכיבו משוואה דיפרנצילית המתארת את התאוצה כתלות בזמן, פיתרו אותה באופן נומרי שרטטו גרף של מהירות הגוף כתלות בזמן

9. נתון מכל מים גלילי שקוטר בסיסו 20cm. גובה המים במכל הוא 60cm. מנקבים בתחתית המכל חריר בקוטר 7mm. לתוך המכל זורמים מים בקצב קבוע של 400cm3s400\frac{cm^3}{s}.
א. הרכיבו משוואה דיפרנציאלית המתארת את המהירות בה יורדים או עולים פני המים במכל.
ב. לפי המשוואה שהרכבתם על איזה גובה מתיצבים פני המים במכל?
ג. עבור נתוני השאלה, שרטטו גרף של גובה המים במכל כתלות בזמן[removed]
ד. לפי הגרף, על איזה גובה מתיצבים פני המים במכל? השוו לתשובה בסעיף ב.


10. תאוצתו של אצן תחרותי מקיימת את המשוואה הדיפרנציאלית [removed] a=A−vb a = A - \frac v b כאשר [removed] A=12.2msA= 12.2\frac ms ו- [removed] b=0.892sb=0.892 s. בעזרת שיטת אוילר שרטטו גרף של מהירות האצן כתלות בזמן בהנחה שהוא מתחיל ממנוחה. כעבור כמה זמן יעבור אצן זה מרחק של 100 מטר?