פתרון משוואות דיפר’ בעזרת ODE על ידי מטלב

Posted on April 22, 2010

0


נניח שיש לנו משוואה דיפר’ מכל מעלה ואנחנו רוצים לדעת את הערכים של המשוואה במספר מוגדר של נקודות. מטלב נותנת לנו את האופציה לחשב בצורה מיידית את הערכים של המשוואה על ידי שימוש בפקודת ode.

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

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

תאור הפונקציה

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

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

נניח שברשותנו המשוואות:

הפונקציה תראה כך:

function xtag = eqt(t,x)
% here
% x(1) = w
% x(2) = q

xtag(1) = -sin(x(2))
% w’ = sin(q)

xtag(2) = x(1);
% q’ = w
xtag=[xtag(1);xtag(2)]

נקרא לקובץ על שם הפונקציה, כלומר eqt.

הרעיון שעומד מאחורי הפונקציה הוא שאנו ניגשים לפתור את הבעיה הזו כבעיית אלגברה לינארית. אנו מקבלים מידע  מהמשתמש לגבי נקודות הבדיקה הרצויות (x,t), אלה יוצבו בתוך המשוואות הדיפר’ ((xtag(1), xtag(2) ויפלטו חזרה כמטריצה (xtag).

תאור קובץ ההפעלה

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

t0 = 0;

% Start time
tf = 20;
% Stop time
x0 = [0 pi/4]
% Initial conditions
[t,s] = ode45(@eqt,[t0,tf],x0);
x = s(:,1)
y = s(:,2)
plot(t,x,’g’,t,y,’r’)

ניתן לראות שהגדרנו תנאי התחלה חדשים (x0) בצורת מטריצה יחסית לסידור שעשינו בשורה האחרונה בפונקציה שכתבנו בסעיף קודם. כמו כן הוגדרו זמן התחלה וסוף (t0 ו tf) – ניתן למימוש גם על ידי פקודת tspan. כמו כן אנו דואגים לריצת ode על פי פרמטרים שהוכנסו לפונקציה (eqt, t0,tf,x0).

לצורך שירטוט הגרף נבקש לבצע הפרדה בין המידע שקיבלנו חזרה (t,s) מה ode ונבקש plot . את כל זה ניתן לראות בשלושת השורות האחרונות של הקובץ.

כמה נקודות על אוקטב (החלופה החינמית של מטלב)

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

בכל אופן אני מצרף קישור ישיר למידע על הפקודות האלה באוקטב ומאחל לכולנו סופ”ש נעים.

וקסמן

Tagged: