matlab Trajectory of a cannon ball

reistlin9

New member
matlab Trajectory of a cannon ball../images/Emo53.gif

אם מישהו כאן מתמצא במטלאב ? בעיית מסלול פגז תותח ציר XY עם חיכוך הבעיה מוגדרת כאן: http://www8.cs.umu.se/kurser/5DV035/VT08/Matlab/TFyETLab2.pdf הקוד(המבוסס על דפי ההסבר) לכאורה נכון, השתמשתי בODE45 נראה כי התאוצות מחושבות נכון אך תוצאות המרחק המחושב אינן נכונות הפונקציה:
function [dx] = canonball23(t,x,mu) dx=zeros(4,1); g=9.8; dx(1) =x(1); dx(2) =- mu*x(1)*(sqrt(x(1)^2+x(3)^2)); dx(3) =x(3); dx(4) =-g -mu*x(3)*(sqrt(x(1)^2+x(3)^2)); return​
מיו מקדם חיכוך x1 x3 מהירויות ציר X וציר Y dx1 נגזרת ראשונה - מהירות ציר X dx2 נגזרת שניה תאוצה ציר X dx3 נגזרת ראשונה מהירות צירY dx4 נגזרת שניה תאוצה ציר Y אם למישהו יש סקריפט שעובד אשמח לקבל
 

ארול01

New member
יש לך טעות במשוואות ובהגדרות

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

reistlin9

New member
הבעיה שניתבקשנו לפתור אינה זהה

במדויק לבעיה בדף האינטרנט אצלינו מיו קבוע ולכן נכנס כקלט לפונקציה. וכן הוגדר לנו להשתמש ב ODE45 שהוא הסולבר בחירת מחדל - Runge-Kutta ולא אוילר - נחשב טוב יותר משיטת אוילר לגבי הטעות כביכול - זו לא טעות --ODE אינה פותרת משוואות מדרגה שניה ולכן מגדירים משתנה ביניים וכאן הדוגמה בדף האינטרנט נכונה ראה בתמונה מצורפת ה די אקס 1 הינו מהירות (נגזרת ראשונה של מקום) והיא שווה למשתנה הביניים אקס 1 - מהירות בציר אקס האקס 1 מופיע בצד ימין של משוואת הנגזרת השניה - המוצגת כנגזרת ראשונה של המהירות (יתכן כמובן שלא הבנתי אותך ואם זה כך מתנצל מראש- נסה לתקן את הקוד שלי למען ההבנה )
 

ארול01

New member
לא הבנת

function [dx] = canonball23(t,x,mu)zz dx=zeros(4,1);zz g=9.8; dx(1) =x(2);zz dx(2) =- mu*x(2)*(sqrt(x(2)^2+x(4)^2)zz); dx(3) =x(4)zz; dx(4) =-g -mu*x(4)*(sqrt(x(2)^2+x(4)^2)zz); לדעתי זה הקוד הנכון, הוקטור X הוא : x1 מיקום בציר x . x2 מהירות בציר x ולכן שווה גם לנגזרת של 1. x3 מיקום בציר y . x4 מהירות בציר y. אתה עשית סלט מהמשתנים :)
 

reistlin9

New member
המשוואות אינן מכילות תלות מפורשת

בווקטורי המיקום כאן המשוואות הנתונות מסדר שני - נגזרת שניה- התאוצה - תלויה בנגזרת ראשונה המהירות אין תלות מפורשת של משתנה המיקום בנגזרות - ויש קשר בין ציר Y ל X מאחר והתנגדות האוויר מוכפלת בנורמה של ווקטור המהירות. לדוגמה קוד של משוואות לורנץ (מושך לורנץ )שעובד אך כאן המשוואות מראש נתונות ככקשר בין מהירויות לווקטור המיקום http://en.wikipedia.org/wiki/Lorenz_attractor
function [dx] = lorenzo(t,x,ro) dx = zeros(3,1); beta=8/3; sigma=10; dx(1) = sigma*(x(2) - x(1)); dx(2) = x(1)*(ro- x(3))-x(2); dx(3) = x(1)*x(2)-beta*x(3); return​
 

ארול01

New member
שים לב למשוואות שכתבתי

הן העתק מדוייק של המשוואות ברפרנס. 1. הנגזרת הראשונה של המיקום ב x שווה למהירות ב x . 2. נגזרת המהירות ב x שווה לפונקציה של המהירויות ב x ו y . 3. הנגזרת הראשונה של המיקום ב y שווה למהירות ב y . 4. נגזרת המהירות ב y שווה לפונקציה של המהירויות ב x ו y . תנסה, זה עובד.
 

reistlin9

New member
וזה בדיוק מה שכתוב בקוד שכתבתי

DX1 נגזרת של המקום DX2 נגזרת של המהירות DX1 = X1 ו X1 המהירות בציר X dx2=f(x1,x3 המהירויות בציר X ובצירY בהתאמה הפונקציה קולטת את תנאי התחלה דרך משתנה X ומוציאה את הפלט דרך משתנה DX הסולבר משתמש ב משתנה T ומבצע לופ על הפונקציה - בהתאם למהלך האינטגרציה אני לא מצליח לגלות טעות או הבדל בין מה שאתה מציע לשיטה שלי אם ניסית וזה עובד שים כאן את הסקריפט של הפונקציה כולל וקטור תנאי התחלה
 

ארול01

New member
לא הבנת את הרעיון

DX ו X הם משתנים תלויים, הם קשורים בכך ש DX1 זה נגזרת של X1 ולכן אם אתה מגדיר את X1 כמקום ב X אז DX1 זה המהירות והיא לא יכולה להיות שווה ל X1 אלא ל X2 אם אתה מגדיר את X2 כמהירות וכך הלאה בכל המשוואות. אתה צריך שתהיה התאמה בין הוקטור X ו DX ולא להחליט ארביטררית מה כל איבר בכל וקטור. אתה יכול פשוט להעתיק את הפונ' שלי ולהכניס איזה תנאי התחלה שמתחשק לך.
 

reistlin9

New member
את הרעיון הבנתי

MATLAB אינך יודע כנראה משום שאם היית יודע היית שם סקריפט ונגמר העניין
וגם אינך מנסה להבין מה כתבתי אין כאן שום קשר להבט הפיזיקלי של הקשר בן מיקום X לנגזרת של המיקום dx/dt X ו DX הם שמות משתנים - X מייצג מהירות לא מיקום dx מייצג dx/dt נגזרת של המיקום ועל כן גם הוא מהירות. המשוואות עבור התאוצות נכונות ומאטלאב פותרת אותן נכון MATLAB לא מאפשרת לפתור בעיות דיפרנציאליות ממעלה שניה - ודרך הכתיבה הזו היא התחכמות כדי להציג משוואות ממעלה שניה כמשוואה במעלה ראשונה שאותה מאטלאב מסוגלת לפתור. מה לעשות המשוואה הפיזיקלית בנוייה על קשר בן הנזרת השניה לראשונה ומאחר והמשוואות מצומדות ולא לנאריות נח לפתור את הבעיה בצורה נומרית. יש בעיה דומה בציר אחד -בעיית הצנחן גם כאן יש חיכוך והצנחן מגיע למהירות סופית -נפילת טיפות מים - גשם נפתרת באופן דומה - כאן יש תלות בנפח - משקל ושטח החתך של הטיפה הנופלת- עבור הקוד של בעיית הצנחן יש פתרון נכון במאטלאב דוגמה (אם בא לך לראות) יש כאן http://www.youtube.com/watch?v=fx3bl4oA_0U ראה איך בנו את הוקטור RK המקביל לוקטור DX אצלי משום מה הקוד של בעיית פגז תותח - משוואות מצומדות לא עובד. מאטלאב מבצעת נכון את האינטגרציה על התאוצה - משתני DX2 ו DX4 - ומתקבלות מהירויות נכונות - אך את משתנה המיקום מאטלאב אינה פותרת נכון ויש לי כנראה שגיאה היכן שהוא בהצבה למשתני DX1 ו DX3 בקיצר אם מישהו כאן מכיר את הבעיה ויש לו סקריפט של הפונקציה - אשמח לקבלה אם לא לא נורא תודה על העזרה עד כאן יוני
 

ארול01

New member
נתתי לך את הסקריפט

זה שאתה מתעצל לקחת לבדוק שהוא עובד זה בעיה שלך. אתה זה שלא מצליח לכתוב סקריפט שעובד ולא אני. אתה לא צריך ללמד אותי איך פותרים מד"ר במטלב, אני עושה את זה על בסיס יומיומי בעבודה, אם היית מנסה להבין מה כתבתי היית רואה שהסברתי לך בדיוק איפה הטעות שלך. האינדקסים של וקטור X צריכים להיות מתאימים לאינדקסים של וקטור DX כך ש DX1 זה נגזרת של X1 וכו.. .
 

reistlin9

New member
הסקריפט שלי עובד יופי

בלי כל התוספות של ה ZZ ושאר דברים שאינם נחוצים הפונקציה עובדת הבעיה היתה עם העברה לא נכונה של תנאי התחלה - העברתי נגזרת ראשונה ושניה במקום להעביר מקום ונגזרת ראשונה =מהירות תחילית הגרף עבור מהירות 100 מטר שניה זוית 60 ומיו 0.001 נראה בציור המצורף אני עדיין צריך להכניס תנאי עצירה לסולבר שיעצור ולשפר את הפלוטינג בכל מקרה תודה עלהנכונות לעזור
 
למעלה