assignment13_numerical solution of differential equations.ipynb

Kernel: Python 3 (Anaconda 2020)

import matplotlib.pyplot as plt
import numpy as np

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

Shlomo Rozenfeld

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

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

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

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

קטע הקוד שבהמשך מתארת את החישוב עבור משוואה דיפרנציאלית מהצורה [removed]$\frac{dy}{dx}=x^2$. ידוע כי עבור x=3 ערכו של y הוא 1.0 . הבחירה של [removed]$\Delta x=h$ תלויה בדיוק הרצוי. ככל שנבחר $\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]$\frac{dy}{dx}=x^2$ עבור תנאי ההתחלה y(3)=1 פיתרון אנאליטי: [removed]$y(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]$\frac{dN}{dt}=-5N$. כמו כן ידוע כי עבור זמן t=0, N(0)=6.0.

א. שרטטו גרף של פתרון נומרי של המשוואה עבור? [removed]$0\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. מה המסקנה?

הערה: השגיאה היחסית מוגדרת כ- $\frac{|N_{numeric}-N_{true}|}{N_{true}}$

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

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

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

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

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

תרגילים

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

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. לתוך המכל זורמים מים בקצב קבוע של $400\frac{cm^3}{s}$.
א. הרכיבו משוואה דיפרנציאלית המתארת את המהירות בה יורדים או עולים פני המים במכל.
ב. לפי המשוואה שהרכבתם על איזה גובה מתיצבים פני המים במכל?
ג. עבור נתוני השאלה, שרטטו גרף של גובה המים במכל כתלות בזמן[removed]
ד. לפי הגרף, על איזה גובה מתיצבים פני המים במכל? השוו לתשובה בסעיף ב.

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