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

מבוא

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

שיטת אויילר

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

כאשר נתון, ונרצה למצוא את .

פיתוח השיטה

נגדיר את הנקודות:

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

נציב במד”ר, ונקבל ביטוי ל-:

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

כאשר הוא התנאי התחלה שלנו.

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

השגיאה (האמיתית) של שיטה זו וכל שיטה לפתרון נומרי של מד”ר נקראת השגיאה הכללית/גלובלית, והיא נתונה ע”י:

נגדיר עוד שגיאה לשיטות נומריות לפתרון מד”ר:

הגדרה:

שגיאה מקומית . שגיאה מקומית היא השגיאה שנוצרת בצעד יחיד של השיטה. כלומר:

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

כאשר . נסיק כי לשיטה זו שגיאה מקומית מסדר .

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

כאשר בקטע, ו- בקטע.

דוגמה:

נביט במד”ר הבא:

הפתרון האמיתי הוא . נפתור בשיטת אויילר:

נקבל את הפתרונות הבאים עבור ( הוא הפתרון האמיתי, ואילו הוא הפתרון המשוער):

עבור :

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

שיטות רונגה-קוטה (Runge-Kutta)

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

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

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

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

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

שיטת אויילר המתוקנת

אם נבצע אינטגרציה על המד”ר שלנו לפי , בקטע , נקבל:

לפי המשפט היסודי של החדוא, נקבל:

נעביר אגפים:

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

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

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

נשים לב כי שיטה זו כוללת שני חישובים של בשלב . נאמר כי זוהי שיטת בשני שלבים.

השגיאה המקומית:

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

שיטת אויילר המשופרת

נוכל גם לקבל קירוב עבור האינטגרל בעזרת קירוב ביניים, כדי לקבל את שיטת הביניים הסתומה:

כאשר:

נשתמש בשיטת אויילר כדי למצוא קירוב ל-:

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

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

שיטת רונגה-קוטה הקלאסית

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

הפיתוח של השיטה וההוכחה שהשגיאה של היא מסדר היא מסובכת שאפילו לספר לא היה כוח להראות אותו.

נהוג לסמן , ואז השיטה הופכת להיות:

להלן פונקציית MATLAB הכתובה עבור מערכת מד”ר. שימו לב שבמקום לחשב את , מחושב הערך :

function [t,y] = rk4(f,tspan,y0,h)
	%
	% function [t,y] = rk4(f,tspan,y0,h)
	%
	% A simple integration routine to solve the
	% initial value ODE y’ = f(t,y), y(a) = y0,
	% using the classical 4-stage Runge-Kutta method
	% with a fixed step size h.
	% tspan = [a b] is the integration interval.
	% Note that y and f can be vector functions
	
	y0 = y0(:); % make sure y0 is a column vector
	m = length(y0); % problem size
	t = tspan(1):h:tspan(2); % output abscissae
	N = length(t)-1; % number of steps
	y = zeros(m,N+1);
	y(:,1) = y0; % initialize
	
	% Integrate
	for i=1:N
		% Calculate the four stages
		K1 = feval(f, t(i),y(:,i) );
		K2 = feval(f, t(i)+.5*h, y(:,i)+.5*h*K1);
		K3 = feval(f, t(i)+.5*h, y(:,i)+.5*h*K2);
		K4 = feval(f, t(i)+h, y(:,i)+h*K3 );
	
	% Evaluate approximate solution at next step
		y(:,i+1) = y(:,i) + h/6 *(K1+2*K2+2*K3+K4);
	end
end

דוגמה:

נביט במד”ר הבא:

הפתרון האמיתי הוא . נראה את השגיאות ב-:

ניתן באמת לראות שכל השיטות נותנות שגיאות מסדר בהתאמה.

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

כאשר נתונה לנו מערכת משוואות דיפרנציאליות עם תנאי התחלה:

נוכל לכתוב אותה בצורה הוקטורית הבאה:

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

דוגמה:

נפעיל את שיטת שיטת רונגה-קוטה הקלאסית על הבעיה הבאה:

נמצא את הפתרון בקטע , עם גודל צעד . קוד MATLAB:

y0 = [80,30]’; % initial data
tspan = [0,100]; % integration interval
h = .01; % constant step size
[tout,yout] = rk4(@func,tspan,y0,h);
figure(1)
plot(tout,yout)
xlabel(’t’)
ylabel(’y’)
legend(’y_1’,’y_2’)
 
figure(2)
plot(yout(1,:),yout(2,:))
xlabel(’y_1’)
ylabel(’y_2’)
 
function f = func(t,y)
  a = .25; b = -.01; c = -1; d = .01;
  f(1) = a*y(1) + b*y(1)*y(2);
  f(2) = c*y(2) + d*y(1)*y(2);

הפתרון מוצג באיור הבא:

הפלט של הבעיה הנ”ל.

---

תרגיל:
נתונה המשוואה הדיפרנציאלית הבאה מסדר עם תנאי התחלה:

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

פתרון:
נבצע הורדת הסדר:

נציב במשוואה:

לכן הבעיה:

בצורה וקטורית:

---

מד”ר עם תנאי שפה

בכללי, מערכת מד”ר עם משוואות הנתונה בצורה הבאה

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

דוגמה:

את המד”ר:

נוכל לרשום כמערכת משוואות:

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

אז למשוואה יש כעת שני פתרונות שונים:

שני פתרונות שונים לבעיה עם תנאי שפה.

שיטת הירי (Shooting method)

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

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

אנו נתמקד בפתרון נומרי של בעיה מסדר :

נפרק לשני מקרים - לינארית ו- לא לינארית.

שיטת הירי עבור משוואה לינארית

נביט במקרה הכללי של בעיה לינארית מסדר עם תנאי שפה:

כיוון שהבעיה לינארית, נוכל למצוא פתרון עבורה בעזרת רק “ניחושים” - כלומר רק שני בעיות התחלה () עם תנאי התחלה שונים. פעם אנו “יורים” בשיפוע , פעם בשיפוע :
בעיה :

בעיה :

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

נדרוש שהפתרון שלה יקיים את תנאי השפה הראשון, :

מהגדרת הבעיות, ולכן:

נדרוש שהפתרון שלה יקיים את תנאי השפה השני, :

מהסופרפוזיציה שפיתחנו מקודם:

אין כבר משמעות ל- ב-, אז פשוט נסמן אותו . נסכם:

הערות:

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

---

תרגיל:
פתור באמצעות שיטת הירי את הבעיה הנ”ל:

פתרון:
נרצה למצוא את :

נבחר שיפועים התחלתיים:

ונקבל מהפתרון בשיטת רונגה-קוטה הקלאסית עם ש:

נציב בחזרה בנוסחה עבור ונקבל כי:

ולפיכך:

קוד python:

import numpy as np
import matplotlib.pyplot as plt
 
# define the Runge-Kutta 4th order method for second order ODEs, that returns the solution at each step
def rk4(f, y0, t0, t1, h):
	t = t0
	y = y0
	sol = []
	while t < t1:
		k1 = h * f(t, y)
		k2 = h * f(t + h/2, y + k1/2)
		k3 = h * f(t + h/2, y + k2/2)
		k4 = h * f(t + h, y + k3)
		y = y + (k1 + 2*k2 + 2*k3 + k4)/6
		t = t + h
		sol.append(y)
	return np.array(sol)
 
# define the system of ODEs y''+y=0
def f(t, y):
	return np.array([y[1], -y[0]])
 
# Initial guesses for y'(0) = -1 and y'(0) = 2
y1_slope = -1
y2_slope = 2
 
# Define the interval of integration and the step size
t0 = 0
t1 = np.pi/2
h = np.pi/50
 
# Solve the system of ODEs using the Runge-Kutta 4th order method
y1 = rk4(f, np.array([0, y1_slope]), t0, t1, h)
y2 = rk4(f, np.array([0, y2_slope]), t0, t1, h)
 
# Linear superposition of the two solutions to get the solution of the ODE with y'(0) = 1
l = (y2_slope - 1)/(y2_slope - y1_slope)
y = l*y1 + (1-l)*y2
 
# The real solution of the ODE is y(t) = sin(t)
t = np.arange(t0, t1, h)
y_real = np.sin(t)
 
# Plot the y1, y2m the numerical solution and the real solution
plt.plot(t, y1[:, 0], label='y1')
plt.plot(t, y2[:, 0], label='y2')
plt.plot(t, y[:, 0], label='y: Numerical solution')
plt.plot(t, y_real, label='sint: Real solution', linestyle='--')
plt.xlabel('t')
plt.ylabel('y')
plt.title('Solution of y\'\' + y = 0')
plt.legend()

הפלט:

הפתרונות , הפתרון הנומרי והאמיתי לבעיה , עם תנאי השפה והניחושים הנ”ל.

---

שיטת הירי עבור משוואה לא לינארית

נביט במקרה הכללי של בעיה לא לינארית מסדר עם תנאי שפה:

נגדיר את הבעיה , שתפתור את בעיית ההתחלה:

אנו רוצים למצוא עבור איזה שיפוע “ירי” נקבל ש- . כלומר, אנו נידרש למצוא שורש של המשוואה:

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