אימייל

צבעים

צילומי נופים

Virtual editor

מילון זעיר

היווצרות החיים

סילונים אסטרונומיים Fulltext

סופרנובה

טורנדו

יחסות IARD

מאמרים על יחסות

ספר על יחסות

סימולציה בעזרת מחשב

ד"ר נציבי בן-אמוץ

by Netsivi Ben-Amots, Haifa, Israel, 1998 Copyright ©

כל הזכויות שמורות לנציבי בן-אמוץ, חיפה, ישראל.

1. זוית מתבדרת זה בכלל לא בידור

מצאתי פתרון אנליטי לזוית הפרצסיה של רוטור בעל סימטריה סיבובית, ונסיתי לאמת אותו בסימולציה במחשב. השתמשתי בשפת המחשב CSMP, שהיתה נוחה מאוד, והמשוכללת מסוגה באותה תקופה (1970). את צעד זמן האינטגרציה בחרתי כ-עשירית מזמן המחזור של התדירות הגבוהה ביותר, שהיה ידוע לי מהפתרון האנליטי שמצאתי.

כשהרצתי את הסימולציות הראשונות במחשב, הזוית טתה (θ) התבדרה. התיעצתי עם הסטודנט א.מ. שנחשב למומחה. א.מ. הסביר לי שזו בעיה נומרית המופיעה במקרים מסוימים. הפתרון הוא להוסיף פקודה, שאם הזוית עולה על 180°, יש להחליף אותה בזוית שקטנה ממנה מ-180°. עושים זאת ע"י פונקציה הנקראת AMOD, המחשבת את השארית (ה"מודולו") מחלוקת הזוית ב-180. מובן שאם החישוב הוא ברדיאנים ולא במעלות, מחלקים ב-π (ז"א ב-3.1416) במקום ב-180. זה פתר את הבעיה.

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

את הסימולציות במחשב עשיתי בשפת CSMP, שהיתה נוחה במיוחד לתיכנות סימוליציה במחשב. מבין אפשרויות האינטגרציה, בחרתי בזו שהיתה המדויקת ביותר, שיטת Runge-Kutta מסדר רביעי. בשיטה זו תכנת CSMP היתה משנה את מירווח הדגימה (delt) בשיטה מתוחכמת הנקראת predictor-corrector. בשיטה זו התכנית בודקת את תוצאות הסימולציה, ואם הן נראות "חשודות", התכנית מקטינה את מירווח האינטגרציה. מדי פעם התכנית בודקת אם לדעתה אפשר להגדיל את מירווח האינטגרציה delt, וע"י כך לחסוך זמן מחשב. חשדתי בשיטה זו, ורציתי לבדוק אותה. בשיטה זו אם התכנית מצאה שיש להקטין את מירווח האינטגרציה (delt), היא חזרה אחורה בזמן time, והמשיכה משם באינטגרציה עם מירווח אינטגרציה delt מוקטן.

כדי לבדוק, שמרתי את זמן האינטגרציה time במשתנה נפרד. בחישוב הבא בדקתי אם הזמן time גדול או קטן מאשר היה בחישוב הקודם. אם הוא קטן מהזמן הקודם, זה סימן שה-predictor-corrector פעל, והחזיר את הזמן אחורה, כדי לחשב מזמן זה ואילך במירווח זמן (delt) יותר קטן. הוספתי פקודה שאם כך קרה, תודפס המילהAGAIN והזמן time בו זה קרה והזמן ש"חזר אחורה". השויתי את הזמנים לתדירות החדשה של התנודות שהמחשב מצא, ובאמת זו היתה אותה תדירות. התנודות היו נומריות, ללא קשר עם המשוואות.

הפתרון היה פשוט: במקום להפעיל את אינטגרצית Runge-Kutta עם ה-predictor-corrector, שזו היתה ברירת המחדל של שפת CSMP, בחרתי באינטגרציה בשיטת Runge-Kutta בצעד קבוע ((RKSFX בשפת CSMP. התיעצתי שוב עם הסטודנט המומחה א. מ., והוא אישר שאכן אינטגרציה במירווח זמן משתנה יכולה לגרום לתנודות נומריות.

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

היועץ ד"ר מנחם שמואלי ז"ל (1935-1980), שכבר עשה סימולציות במחשב בשפת פורטרן למחקריו, אמר לי שהדרך היחידה היא להגדיל את דיוק החישוב, ולראות מה יקרה. הגדלת דיוק החישוב (כמה ספרות "אחרי הנקודה העשרונית") היתה כרוכה בעבודה ענקית: בדיוק הרגיל עמדה לרשותי שפת המחשב CSMP, שהיתה נוחה מאוד, אבל באותה תקופה איפשרה רק דיוק רגיל.

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

נאלצתי לכתוב תוכנית סימולציה בשפת FORTRAN עם דיוק כפול.

פניתי שוב לסטודנט א.מ., שנתן לי דוגמה לסברוטינה בשפת פורטרן, המבצעת אינטגרציה בשיטת Runge-Kutta בצעד קבוע בשפת FORTRAN. עבדתי הרבה על התאמת התכנית לצרכי. בניתי את התכנית כך שתהיה דומה ל-CSMP, עם שלושה חלקים: INITIAL .1 בו עושים חישובים מוקדמים, ובו הגדרתי את מירווח האינטגרציה delt וזמן סיום האינטגרציה fintim ואת תנאי ההתחלה וכן הלאה, וחישבתי את תנאי ההתחלה שרציתי לבדוק, 2. DYNAMIC, ובו יכולתי לפרט את המשואות הדינמיות, ו-3. TERMINAL בו עשיתי חישובים סופיים וקריאה לסברוטינות השירטוט. הוספתי את הדיוק הכפול. שלא כמו בשפת CSMP, שמרתי את דגימות תוצאות החישובים בזכרון המחשב, ורק בתום החישוב דגמתי מזכרון המחשב את הנקודות לצורכי הגרפים. כך השגתי חסכון משמעותי בזמן המחשב שנדרש לכל ריצה (ולא עם PREDICTOR-CORRECTOR שמזייף את התוצאות). החסרון בשיטה זו הוא שאם הריצה נכשלת, אין שמירה של התוצאות עד לרגע הכשלון, ויותר קשה למצוא מה הכשיל את ריצת המחשב. אבל זו רק יותר עבודה במציאת השגיאות שלי, אבל לא פגיעה באמינות תוצאות המחשב אם התקבלו. מה לעשות, בשביל תוצאות טובות צריך לעבוד יותר...

בנוסף לגרפים של המשתנים לפי הזמן, ציירתי גם גרף פולרי של הזוית תטה θ לפי הרדיוס, וכך קיבלתי ספירלה מתבדרת או מתכנסת – מה שלא היה אפשרי בשפת CSMP. לשם כך פיתחתי סברוטינה "אוטומטית" לשירטוט משתנה אחד לעומת משתנה שני שאינו זמן. זה לא היה פשוט: CSMP וכל סברוטינה אחרת הדפיסו גרפים במדפסת. גישה ל-PLOTTER היתה קשה, יקרה ולא מעשית לאלפי הריצות שעשיתי. אבל בהדפסה במדפסות שהיו אז, היה כלל ברזל: אחרי שהמדפסת הדפיסה שורה והתקדמה לשורה הבאה, המדפסת לא יכלה לחזור אחורה. אי אפשר היה לצייר ספירלה נקודה אחרי נקודה. לכן, חישבתי מראש מערך דו-מימדי השומר בזכרון המחשב את כל נקודות הגרף. רק כשנגמר חישוב כל הנקודות המערך הדו-מימדי, יכולתי להדפיס ממנו את הגרף במדפסת, שורה אחרי שורה, ולקבל ספירלה מתכנסת או מתבדרת או אליפטית וכו', בהתאם למשואות הדינמיות.

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

עצה חשובה קיבלתי מהמנחה ד"ר יצחק פורת ז"ל (1934-2012). הוא הסביר שתכנית מחשב יש לבדוק תחילה על מקרה\ים בו\בהם יש פתרון אנליטי, ולוודא שתוצאות הסימולציה הנומרית מתאימות לפתרון הידוע. רק לאחר מכן יש טעם להשתמש בתכנית המחשב לפתרון מקרים בהם הפתרון האנליטי אינו ידוע. (באופן כללי, הרחבתי כלל זה גם לתיאוריות חדשות – הן צריכות תחילה להגיע לפתרון נכון במקרים בהם הפתרון ידוע בדרך אחרת).

רציתי לבדוק עוד מקרים, וכמי שנכווה ידעתי שכל שאטפל במקרים יותר מסובכים, שאין לי עבורם הערכה אנליטית, וידוע לי פחות, או לא ידוע לאיזה תוצאות עלי לצפות בהם, אני עלול להגיע לאי יציבות נומרית גם בדיוק כפול. כדי לוודא שלא נכשלתי בכך, בכל מקרה חדש הרצתי בהתחלה אינטגרציות בשתי השיטות המדויקות Runge-Kuttaמסדר רביעי ו-Adams מסדר שני, ווידאתי שההפרש בין תוצאות הריצות זניח, ורק אז הרצתי בשיטה המדויקת יותר: Runge-Kutta מסדר רביעי.

מצויד בטכניקות האלו, בדקתי גם מקרים מורכבים יותר:

1. מקרים בהם התלות במשואות לא איפשרה "התקדמות" בזמן כי משתנה היה תלוי ב"עצמו" כך שלא ניתן לחילוץ. CSMP טיפל במקרים כאלו ע"י שאיפשר להגדיר משתנה שאיננו ניתן לחילוץ כ-IMPLICIT, ואז CSMP דאגה לחשב איטרציות בכל צעד זמן בלי להגדיל את הזמן, עד שהתקבלה התכנסות וניתן היה לדגום ולהשתמש במשתנה הסורר.

כשתיכנתתי מקרים כאלו בשפת FORTRAN, הוספתי LOOPS שביצעו פעולה דומה. למעשה "דרכתי במקום" מבחינת הזמן, כאשר ההתקדמות וההתכנסות היו במשתנה אחר, שמבחינה לוגית (בזכרון של המוח בראש שלי), חשבתי ש"כיוון" ההתקדמות הוא בניצב לזמן, על גרף דמיוני דו-מימדי.

2. במקרים מסוימים, לא היה די בכך. הוספתי LOOPS בהם משתנה הזמן נשאר קבוע, והתקדמתי אחורה וקדימה בזמן, כדי להגיע להתכנסות הנדרשת. מבחינת המחשב לא היה הבדל אם אני מתקדם ב- LOOPS "ניצב לזמן", או על משתנה זמן "נוסף", "קדימה ואחורה" על ציר הזמן, כל עוד שמרתי שהתקדמות האינטגרציה בין ה-LOOPS היתה בצעד קבוע.

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

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

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

תוכניות המחשב כשלעצמן, אינן הפתרון הסופי לכל בעיה.

כמו שצריך להבין את הפיזיקה ו\או המכניקה של הבעיה, יש להבין גם את הבעיות הנומריות שיכולות להופיע בסימולציה נומרית ובאינטגרציה נומרית. אסור לסמוך על התוצאות בעיניים עצומות, ולהסתפק בהסוואת ה"כביסה המלוכלכת". יש להבליט את הכביסה המלוכלכת, ע"י טריקים כמו ה-AGAIN לעיל, ולהתגבר על הבעיות ולפתור אותן, במקום להסוותן.

2. עוד על סימולציה במחשב

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

אם לסדר את שיטות האינטגרציה הנפוצות למשואות דיפרנציאליות רגילות, לפי דיוקן ויציבותן, הסדר מהשיטות הפחות טובות אל היותר טובות, הוא:

א. אוילר

ב. טרפז

ג. סימפסון

ד. אדמס מסדר שני

ה. רונגה קוטה מסדר רביעי

שיטת אינטגרציה בצעד זמן משתנה, כאשר צעד הזמן נקבע בשיטת PREDICTOR-CORRECTOR נחשבת מודרנית, אבל ראינו לעיל שאינטגרציה בצעד זמן משתנה, כאשר הצעד נקבע בשיטת PREDICTOR-CORRECTOR, יכולה להביא לתנודות נומריות שאין להן קשר למציאות.

מה אני חושב כאשר אני קורא על הסימולציות של תיאורית ההתנגשות שבה לכאורה נוצר הירח, שאומתה ע"י סימולציות דינמיות רבות במחשב-על, עם מערכות משואות דיפרנציאליות רגילות, בשיטת אינטגרציה של אוילר, כאשר הצעד נקבע בשיטת ?PREDICTOR-CORRECTOR לי זה נראה ביזבוז של מחשב-על, שלא לדבר על פתרון ההתנגשות, בו יציבות הפתרון לא הוכחה ולא נסתרה למרות המאמץ העצום שהושקע כדי להשתמש במחשב-על בסימולציות אלה.

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

אלו כן אלו הן וריאציות מתוחכמות של האימרה הידועה בקשר למחשבים: Garbage in – garbage out

.3. Bio-Nedical Database

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

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

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

עיינתי בהדפסה, ופתאום הבנתי. היא הזדקקה לתכנית הסטטיסטיקה Bio-Medical Database, אבל הדפיסה BND במקום BMD. אין פלא שהמחשב לא הבין מה רוצים ממנו. אבל גם לאדם קשה למצוא טעות כזו, הן מפני שהאותיות M ו-N כל כך דומות, והן מפני שלא מחפשים טעות בשם, ומי שבודק יכול שלא לדעת שלצירוף האותיות BMD יש משמעות, ולכן צירוף אותיות דומה לו הוא הטעות.

מראי מקום

הקורא מוזמן לעיין במראי המקום של עבודות שעשיתי על רוטורים מהירים תוך שימוש בסימולציה במחשב כמתואר לעיל:

1. נ. בן-אמוץ, "תנועת רוטור מהיר בהשפעת המומנט הניצב לצירי הפרצסיה והנוטציה", תיזה ל-.Dsc (דוקטורט), הטכניון, חיפה, ישראל (1976).
2. N. Ben-Amots,

"Approximate analytical solution for high-speed spin-axisymmetric rotor, using coordinate system linked to precession and nutation"

Acta Mechanica, v. 25, No. 1-2, pp. 111-119 (March 1976).
https://doi.org/10.1007/BF01176934

בחזרה לדף הקודם

מעבר לעבודות אחרות:

1. הוראה

2. תכנון צנרת

עבודות או ארועים שכללו תוכנת מחשב או שימושים במחשב (11 – 3):

3. כמעט הייתי HACKER

4. סימולציה נומרית במחשב

5. תכניות מחשב במנהרות הרוח

6. המחשב לעזרת מגדלי הזיקוק

7. המחשב המובטל

8. זמן אמיתי

9. שיפורים במחשבי הנדסה ביו-רפואית

10. יומיים כמהנדס תוכנה בחברה פרטית  (מתוך שנה וחצי) 

11. מומחה לשעבר

12. הכנסות של מיליארדים

חזרה לדף הקודם

קישורים לחלק מהאפשרויות הנוספות באתר זה Relativity, gravitation and relativistic rotation Clarifying some paradoxes of relativity at the extreme ספר על יחסות ---------> להזמנת הספר שלח אימייל אל dillone.bookorder@gmail.com טורנדו מכניזם הטורנדו ו-מקור האנרגיה של הטורנדו מקור החיים על כדור הארץ גלי ים נשברים סופרנובה: כוכב מתפוצץ איך קריסת כוכב הופכת להתפוצצותו סילונים אסטרונומיים  שלושה מאמרים על סיבוב יחסותי:   שלושה מאמרים נוספים על יחסות:  רשימת מאמרים מאת נציבי בן-אמוץ מילון זעיר לתייר דובר אנגלית:  איטליה יוון צרפת ישראל צילומי נופים כתבן וירטואלי כתבן בלי שמירה בדיסק אבל זמין באינטרנט הכנסה של מיליארדים למדינת ישראל IARD International Association for Relativistic Dynamics בחירת צבעים

להזמנת הספר שלח אימייל אל dillone.bookorder@gmail.com הניע את העכבר מעל התמונה. לפרטים נוספים לחץ עם העכבר השמאלי לחץ למטה עם העכבר השמאלי עבור דף הספר באנגלית אצל המו"ל בחו"ל