روش رانگ کوتای مرتبه دوم — از صفر تا صد
در آموزشهای پیشین مجله فرادرس، با روش رانگ کوتای مرتبه اول و مرتبه چهارم آشنا شدیم. در این آموزش، روش رانگ کوتای مرتبه دوم (Second Order Runge-Kutta) را برای حل معادلات دیفرانسیل معرفی خواهیم کرد.
تعریف مسئله
معادله دیفرانسیل زیر با جواب $$ y ( t) $$ را در نظر بگیرید:
$$ \large { { d y ( t ) } \over { d t } } = y' ( t ) = f ( y ( t ) , t ) , \quad \quad { \rm{\;}} y ( t _0 ) = y _ 0 $$
میخواهیم جواب $$ y ( t) $$ را با کامپیوتر تقریب بزنیم (با دانستن شرایط اولیه $$ y ( t _ 0 ) = y _ 0 $$). در ادامه، یک روش شهودی برای انجام این کار ارائه کرده و مثالهایی را نیز بررسی میکنیم. این روش همان روش رانگ کوتای مرتبه دوم است.
روش رانگ کوتای مرتبه دوم
در این بخش، روش رانگ کوتای مرتبه دوم را برای دستههایی از معادلات دیفرانسیل به صورت شهودی معرفی خواهیم کرد.
معادله دیفرانسیل مرتبه اول بدون ورودی
همانطور که در آموزشهای پیشین گفتیم، در روش رانگ کوتای مرتبه اول، از مشتق در زمان $$ t _ 0 $$ (در شکل زیر $$ t _ 0 = 0 $$) برای تخمین مقدار تابع در یک گام زمانی بعد استفاده میشود. برای مثال، معادله دیفرانسیل زیر را در نظر بگیرید:
$$ \large { { d y ( t ) } \over { d t } } + 2 y( t ) = 0 $$ یا $$ \large { { d y ( t ) } \over { d t } } = - 2 y ( t ) $$
که شرایط اولیه آن $$ y ( 0 ) = 3 $$ است. جواب دقیق این معادله $$ y ( t) = 3 e ^ {-2t} $$ برای $$ t \ge 0 $$ است که فرض میکنیم آن را نمیدانیم و برای ارائه یک تقریب به انتگرالگیری عددی نیاز داریم.
در شکل زیر، شیب در $$ t = 0 $$ را $$ k _ 1 $$ و تخمین را $$ y^* ( h ) $$ مینامیم که در این مثال، $$ h = 0.2 $$ است.
بدیهی است که مقداری خطای تخمین وجود دارد که میخواهیم آن را کم کنیم. یکی از راههای نیل به این هدف، استفاده از مشتق در نقطه میانی بین $$ t = 0 $$ و $$ t = h = 0.2 $$ است. شیب در این نقطه ($$ t = \frac{1}{2}h=0.1$$) در شکل زیر نشان داده شده و با $$ k _ 2 $$ مشخص شده است. خط نارنجی نیز، در واقع خط مماس بر منحنی آبی در $$ t = \frac{1}{2} h $$ است.
اکنون اگر از شیب میانی $$ k _ 2 $$ استفاده کنیم، تخمین بهتری نسبت به قبل برای $$ y ^ * ( h ) $$ خواهیم داشت. در شکل زیر، مقدار دقیق جواب $$ y (0.1) = 2.0110 $$ و همچنین، تقریب $$ y ^ * (0.1) = 2.0175 $$ برای خطایی در حدود $$ 0.3$$ درصد نشان داده شدهاند. دقت کنید که در روش رانگ کوتای مرتبه اول، این خطا تقریباً $$ 10 $$ درصد است.
به نظر میرسد جواب تقریب بسیار مناسبی دارد و بدیهی است که دقت آن نسبت به جواب به دست آمده از روش رانگ کوتای مرتبه اول (که در آن از خطی با شیب $$ k _ 1 $$ در $$ t = 0 $$ استفاده میشد) بسیار بیشتر است. مسئلهای که با آن مواجهیم، این است که مقدار دقیق $$ y ( \frac{1}{2} h ) $$ را نمیدانیم و به همین دلیل نمیتوانیم مقدار دقیق $$ k _ 2 $$ را بیابیم. لازم به یادآوری است که محاسبه مشتق نیازمند دانستن مقدار تابع است ($$y' ( t ) = - 2 y ( t) $$).
یک راه برای حل این مشکل، استفاده از روش اویلر یا همان روش رانگ کوتای مرتبه اول برای به دست آوردن تقریب $$ y ( t) $$ در $$ t = \frac{1}{2} h = 0.1 $$ است که آن را $$ y _ 1 ( \frac{1}{2} h ) $$ مینامیم. سپس میتوانیم از این تخمین برای به دست آوردن $$ k _ 2 $$ (که تقریبی از شیب در نقطه میانی است) استفاده کنیم و در نتیجه، از $$ k _ 2 $$ برای یافتن $$ y ^* ( t) $$ بهره ببریم. برای رفتن از نقطه شروع در $$ t = 0 $$ تا تخمین در $$ t = h $$، فرایند زیر را انجام میدهیم:
عبارت مشتق در $$ t = 0 $$: $$ \large y' ( 0 ) = - 2 y ( 0) $$
مشتق در $$ t = 0 $$: $$ \large { k _ 1 } = - 2 y ( 0 ) $$
تخمین میانی تابع در $$ t = h / 2 $$: $$ \large { y _ 1 } \left ( { { h \over 2 } } \right ) = y ( 0 ) + { k _ 1 } { h \over 2 } $$
تخمین شیب در $$ t = h / 2 $$: $$ \large { k _ 2 } = - 2 { y _ 1 } \left ( { { h \over 2 } } \right ) $$
سری تیلور حول $$ t = 0 $$: $$ \large y ( h ) = y ( 0 ) + y' ( 0 ) h + y^{\prime \prime } ( 0 ) { { { h ^ 2 } } \over 2 } + \cdots $$
سری تیلور بریده: $$ \large y ( t ) \approx y ( 0 ) + y' ( 0 ) h $$
تخمین $$ y ( h) $$: $$ \large { y ^ * } ( h ) = y ( 0 ) + { k _ 2 } h $$
در حالت کلی، یک گام به جلو از تخمین $$ t = t _ 0 $$ تا $$ t = t _ 0 + h $$ به صورت زیر انجام میشود:
عبارت مشتق در $$ t = t _ 0 $$: $$ \large y' ( t _ 0 ) = - 2 y ( t _0) $$
مشتق در $$ t = t _ 0 $$: $$ \large { k _ 1 } = - 2 y ( t _ 0 ) $$
تخمین میانی تابع در $$ t = t _ 0 + h / 2 $$: $$ \large { y _ 1 } \left ( { {t _ 0 + h / 2 } } \right ) = y ( t_0 ) + { k _ 1 } { h \over 2 } $$
تخمین شیب در $$ t = t _ 0 + h / 2 $$: $$ \large { k _ 2 } = - 2 { y _ 1 } \left ( { { t _ 0 + h / 2 } } \right ) $$
سری تیلور حول $$ t = t _ 0 $$: $$ \large y ( t _ 0 + h ) = y ( t _ 0 ) + y' ( t _ 0 ) h + y^{\prime \prime } ( t _ 0 ) { { { h ^ 2 } } \over 2 } + \cdots $$
سری تیلور بریده: $$ \large y ( t ) \approx y ( t _ 0 ) + y' ( t _ 0 ) h $$
تخمین $$ y ( t _ 0 + h) $$: $$ \large { y ^ * } ( t _ 0 + h ) = y ( t _ 0 ) + { k _ 2 } h $$
تقریب معادله دیفرانسیل مرتبه اول بدون ورودی با استفاده از متلب
میتوانیم از متلب استفاده کنیم و محاسبات بالا را با استفاده از آن پیادهسازی کنیم. روندی که در بالا آن را توضیح دادیم، به صورت زیر در متلب پیاده میشود:
1%% Example 1
2% Solve y'(t)=-2y(t) with y0=3, midpoint method
3y0 = 3; % Initial Condition
4h=0.2; % Time step
5t = 0:h:2; % t goes from 0 to 2 seconds.
6yexact = 3*exp(-2*t) % Exact solution (in general we won't know this
7ystar = zeros(size(t)); % Preallocate array (good coding practice)
8
9ystar(1) = y0; % Initial condition gives solution at t=0.
10for i=1:(length(t)-1)
11 k1 = -2*ystar(i) % Approx for y gives approx for deriv
12 y1 = ystar(i)+k1*h/2; % Intermediate value
13 k2 = -2*y1 % Approx deriv at intermediate value.
14 ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y
15end
16plot(t,yexact,t,ystar);
17legend('Exact','Approximate');
هر خط از کد بالا به خوبی مراحل الگوریتم را نشان میدهد. با تغییر اندکی در برنامه، میتوان اثر اندازه $$ h$$ را بر دقت جواب مشاهده کرد. شکل زیر این موضوع را به خوبی نشان میدهد. همانطور که میبینیم، هرچه مقدار $$ h$$ بزرگتر باشد، تقریب ضعیفتر خواهد بود.
الگوریتم روش رانگ کوتای مرتبه دوم (نقطه میانی)
معادله دیفرانسیل مرتبه اول زیر داده شده است:
$$ \large { { d y ( t ) } \over { d t } } = f ( y ( t ) , t ) $$
برای حرکت رو به جلوی یک نقطه در $$ t = t _ 0$$، یعنی $$ y ^ * ( t _ 0 ) $$، به اندازه یک گام زمانی $$ h$$، مراحل زیر را طی میکنیم:
تخمین مشتق در $$ t = t _ 0 $$: $$ \large { k _ 1 } = f ( { y ^ * } ( { t _ 0 } ) , \; { t _ 0 } ) $$
تخمین میانی تابع در $$ t = t _ 0 + \frac { h } { 2 } $$: $$ \large { y _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) = { y ^ * } ( { t _ 0 } ) + { k _ 1 } { h \over 2 } $$
تخمین شیب در $$ t = t _ 0 + \frac {h } { 2 } $$: $$ \large { k _ 2 } = f \left ( { { y _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) , \; { t _ 0 } + { h \over 2 } } \right ) $$
تخمین $$ y ( t _ 0 + h ) $$: $$ \large { y ^ * } \left ( { { t _ 0 } + h } \right ) = { y ^ * } ( { t _ 0 } ) + { k _ 2 } h $$
دقت کنید که مقدار اولیه یا همان شرایط اولیه در ابتدای الگوریتم مشخص باشد.
از الگوریتم رانگ کوتای مرتبه دوم اغلب با نام الگوریتم نقطه میانی (Midpoint Algorithm) یاد میشود، زیرا از شیب $$ k_ 2 $$ در نقطه میانی استفاده میکند.
معادله دیفرانسیل مرتبه اول با ورودی
با اضافه شدن تابع ورودی به معادله دیفرانسیل، مشکل چندانی به وجود نخواهد آمد. فقط لازم است مقدار تابع در $$ t = t _ 0 $$ را برای یافتن $$ k _ 1 $$ و همچنین، مقدار آن را در $$ t = t _ 0 + \frac { 1 } { 2 } h $$ برای یافتن $$ k _ 2 $$ تضمین کنیم. برای مثال، ورودی $$ \cos ( 4 t ) $$ را در نظر بگیرید.
$$ \large y' ( t ) + 2 y ( t ) = \cos ( 4 t ) $$ یا $$ \large y' ( t ) = - 2 y ( t ) + \cos ( 4 t ) $$
بنابراین، مراحل زیر را طی میکنیم:
عبارت مشتق در $$ t = t _ 0 $$: $$ \large y' ( { t _ 0 } ) = - 2 y ( { t _ 0 } ) + \cos ( 4 { t _ 0 } ) $$
تقریب مشتق در $$ t= t _ 0 $$: $$ \large { k _ 1 } = - 2 { y ^ * } ( { t _ 0 } ) + \cos ( 4 { t _ 0 } ) $$
تخمین میانی تابع در $$ t = t _ 0 + h / 2 $$: $$ \large { y _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) = { y ^ * } ( { t _ 0 } ) + { k _ 1 } { h \over 2 } $$
تخمین شیب در $$ t = t _ 0 + h / 2 $$: $$ \large { k _ 2 } = - 2 { y _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) + \cos \left ( { { t _ 0 } + { h \over 2 } } \right ) $$
تخمین $$ y ( t _ 0 + h ) $$: $$ \large { y ^ * } ( { t _ 0 } + h ) = y ( { t _ 0 } ) + { k _ 2 } h $$
گفتنی است که جواب دقیق این معادله $$ y ( t ) = 2 . 9 { e ^ { - 2 { \kern 1pt} t } } + 0 . 1 \cos ( 4 t ) + 0 . 2 \sin ( 4 t ) $$ است.
تقریب معادله دیفرانسیل مرتبه اول با ورودی با استفاده از متلب
میتوانیم از متلب برای محاسبات بالا استفاده کنیم. برای این تقریب جدید، باید محاسبات $$ k _ 1 $$ و $$ k _ 2 $$ را با استفاده از مقدار مناسب متغیر زمان تغییر دهیم. برنامه زیر در متلب مربوط به مثال اخیر است.
1%% Example 2
2% Solve y'(t)=-2y(t)-cos(4t) with y0=3
3y0 = 3; % Initial Condition
4h=0.2; % Time step
5t = 0:h:2; % t goes from 0 to 2 seconds.
6yexact = 0.1*cos(4*t)+0.2*sin(4*t)+2.9*exp(-2*t); % Exact Solution
7ystar = zeros(size(t)); % Preallocate array (good coding practice)
8
9ystar(1) = y0; % Initial condition gives solution at t=0.
10for i=1:(length(t)-1)
11 k1 = -2*ystar(i)+cos(4*t(i)); % Approx for y gives approx for deriv
12 y1 = ystar(i)+k1*h/2; % Intermediate value
13 k2 = -2*y1+cos(4*(t(i)+h/2)); % Approx deriv at intermediate value.
14 ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y
15end
16plot(t,yexact,t,ystar);
17legend('Exact','Approximate');
با کمی تغییر در برنامه بالا، میتوانیم تقریب را برای $$ h$$های مختلف رسم کنیم. همانطور که میبینیم، با کوچک شدن $$ h$$، تقریبها دقیقتر خواهند بود.
معادله دیفرانسیل غیرخطی مرتبه اول
در این مورد نیز با مشکل خاصی مواجه نخواهیم بود. در اینجا نیز باید تضمین کنیم که مقدار تابع در $$ t = t _ 0 $$ را برای یافتن $$ k _ 1 $$ و همچنین، مقدار آن را در $$ t = t _ 0 + \frac { 1 } { 2 } h $$ برای یافتن $$ k _ 2 $$ داریم. برای مثال، سیستم زیر را در نظر بگیرید:
$$ \large { { d y ( t ) } \over { d t } } - y ( t ) ( 1 - 2 t ) = 0 , \quad \quad y ( 0 ) = 1 $$
میدانیم که جواب این معادله به صورت زیر است:
$$ \large y ( t ) = { e ^ { t - { t ^ 2 } } } $$
روند عددی به دست آوردن جواب معادله نیز مطابق مراحل زیر است:
عبارت مشتق در $$ t = t _ 0 $$: $$ \large y' ( { t _ 0 } ) = y ( { t _ 0 } ) ( 1 - 2 { t _ 0 } ) $$
تقریب مشتق در $$ t= t _ 0 $$: $$ \large { k _ 1 } = y ^ * ( { t _ 0 } ) ( 1 - 2 { t _ 0 } ) $$
تخمین میانی تابع در $$ t = t _ 0 + h / 2 $$: $$ \large { y _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) = { y ^ * } ( { t _ 0 } ) + { k _ 1 } { h \over 2 } $$
تخمین شیب در $$ t = t _ 0 + h / 2 $$: $$ \large { k _ 2 } = - 2 { y _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) \left ( { 1 - 2 \left ( { { t _ 0 } + { h \over 2 } } \right ) } \right ) $$
تخمین $$ y ( t _ 0 + h ) $$: $$ \large { y ^ * } ( { t _ 0 } + h ) = y ( { t _ 0 } ) + { k _ 2 } h $$
تقریب معادله دیفرانسیل غیرخطی مرتبه اول با ورودی با استفاده از متلب
مانند مثالهای قبل، به سادگی میتوانیم برنامه تقریب معادله دیفرانسیل غیرخطی مرتبه اول با ورودی را در متلب بنویسیم. این برنامه به صورت زیر است:
1%% Example 3
2% Solve y'(t)=y(t)(1-2t) with y0=1
3y0 = 1; % Initial Condition
4h=0.2; % Time step
5t = 0:h:2; % t goes from 0 to 2 seconds.
6% Exact solution, hard to find in this case (in general we won't have it)
7yexact = exp(t-t.^2);
8ystar = zeros(size(t)); % Preallocate array (good coding practice)
9
10ystar(1) = y0; % Initial condition gives solution at t=0.
11for i=1:(length(t)-1)
12 k1 = ystar(i)*(1-2*t(i)); % Approx for y gives approx for deriv
13 y1 = ystar(i)+k1*h/2; % Intermediate value
14 k2 = y1*(1-2*(t(i)+h/2)); % Approx deriv at intermediate value.
15 ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y
16end
17plot(t,yexact,t,ystar);
18legend('Exact','Approximate');
شکل زیر، جواب دقیق و تقریبی معادله را به ازای $$h$$های مختلف نشان میدهد.
همانطور که مشاهده میکنیم، با کاهش $$ h$$، دقت تقریب زیاد میشود.
معادله دیفرانسیل خطی مرتبه بالاتر
معادله دیفرانسیل مرتبه بالاتر زیر را در نظر بگیرید:
$$ \large \displaylines {
{ { { d ^ 3 } y ( t ) } \over { d t } } + 4 { { { d ^ 2 } y ( t ) } \over { d t } } + 6 { { d y ( t ) } \over { d t } } + 4 y ( t ) = \gamma ( t ) \quad \quad \quad \gamma ( t ) = { \rm{unit step function}} \cr
{ \left . { { { { d ^ 2 } y ( t ) } \over { d t } } } \right | _ { t = { 0 ^ + } } } = 0 , \quad \quad { \left. { { { d y ( t ) } \over { d t } } } \right | _ { t = { 0 ^ + } } } = - 1 , \quad \quad y ({ 0 ^ + } ) = 0 \quad \cr } $$
معادله بالا را میتوانیم به صورت دستگاه معادلات دیفرانسیل ماتریسی بازنویسی کنیم:
$$ \large \eqalign {
{ \bf { q } }' ( t ) & = \left[ { \matrix {
{ { q _ 1 }' ( t ) } \cr
{ { q _ 2 }' ( t ) } \cr
{ { q _ 3 }' ( t ) } \cr } } \right] = \left[ { \matrix {
0 & 1 & 0 \cr
0 & 0 & 1 \cr
{ - 4 } & { - 6 } & { - 4 } \cr } } \right] \left [ { \matrix {
{ { q _ 1 } ( t ) } \cr
{ { q _ 2 } ( t ) } \cr
{ { q _ 3 } ( t ) } \cr } } \right] + \left[ { \matrix {
0 \cr
0 \cr
1 \cr } } \right ] \gamma ( t ) \cr
{ \bf { q } }' (t ) & = { \bf { A q } } ( t ) + { \bf { B \gamma } } ( t ) \cr
{ \bf { A } } & = \left [ { \matrix {
0 & 1 & 0 \cr
0 & 0 & 1 \cr
{ - 4} & { - 6} & { - 4} \cr} } \right],\quad \quad \quad { \bf { B } } = \left [ {\matrix{
0 \cr
0 \cr
1 \cr
} } \right] \cr} $$
روند به دست آوردن جواب در این حالت، به صورت زیر است:
عبارت مشتق در $$ t = t _ 0 $$: $$ \large \bf { q }' ( { t _ 0 } ) = { \bf { A q } } ( { t _ o } ) + { \bf { B } } \gamma ( { t _ 0 } ) $$
تقریب مشتق در $$ t= t _ 0 $$: $$ \large { { \bf { k } } _ 1 } = { \bf { A } } { { \bf { q } } ^ * } ( { t _ o } ) + { \bf { B } } \gamma ( { t _ 0 } ) $$
تخمین میانی تابع در $$ t = t _ 0 + h / 2 $$: $$ \large { { \bf { q } } _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) = { { \bf { q } } ^ * } ( { t _ o } ) + { { \bf { k } } _ 1 } { h \over 2 } $$
تخمین شیب در $$ t = t _ 0 + h / 2 $$: $$ \large { { \bf { k } } _ 2 } = { \bf { A } } { { \bf { q } } _ 1 } \left ( { { t _ 0 } + { h \over 2 } } \right ) + { \bf { B } } \gamma \left ( { { t _ 0 } + { h \over 2 } } \right ) $$
تخمین $$ \bf {q} ( t _ 0 + h ) $$: $$ \large { { \bf { q } } ^ * } \left( { { t _ 0 } + h } \right ) = { { \bf { q } } ^ * } ( { t _ o } ) + { { \bf { k } } _ 2 } h $$
جواب دقیق معادله برای $$ t \ge 0 $$ برابر با $$ y ( t ) = { 1 \over 4 } + { { \rm { e } } ^ { - t } } { \mkern 1mu } \left ( { \cos ( t ) - { 5 \over 2 } \sin ( t ) } \right ) - { 5 \over 4 } { { \rm { e } } ^ { - 2 { \kern 1pt} t } } $$ است.
تقریب معادله دیفرانسیل مرتبه سوم با استفاده از متلب
برای تقریب لازم است متغیرهای مناسب اسکالر را با بردار یا ماتریس تغییر دهیم و ماتریسهای $$ \mathbf { A} $$ و $$ \mathbf { B} $$ را تعریف کنیم. توجه کنید که در این حالت، ماتریس $$ \mathbf { B} $$ را در $$1$$ ضرب میکنیم، زیرا ورودی یک پله واحد است ($$ \gamma ( t) = 1 $$ برای $$ t \ge 0 $$). اگر ورودی یک شکل موج سینوسی بود، باید آن در در $$ \mathbf { B} $$ ضرب میکردیم. اگر ورودی وجود نداشته باشد، لازم نیست ماتریس $$ \mathbf { B} $$ را در نظر بگیریم. کد زیر مربوط مثال اخیر است.
1%% Example 4
2% Solve y'''(t)+4y''(t)+6y'(t)+4y(t)=gamma(t), y''(0)=0, y'(0)=-1, y(0)=0
3q0 = [0; -1; 0]; % Initial Condition (vector)
4h=0.2; % Time step
5t = 0:h:5; % t goes from 0 to 2 seconds.
6A = [0 1 0; 0 0 1; -4 -6 -4]; % A Matrix
7B = [0; 0; 1]; % B Matrix
8
9yexact = 1/4 + exp(-t).*(cos(t) - 5*sin(t)/2) - 5*exp(-2*t)/4; % Exact soln
10qstar = zeros(3,length(t)); % Preallocate array (good coding practice)
11
12qstar(:,1) = q0; % Initial condition gives solution at t=0.
13for i=1:(length(t)-1)
14 k1 = A*qstar(:,i)+B*1; % Approx for y gives approx for deriv
15 q1 = qstar(:,i)+k1*h/2; % Intermediate value
16 k2 = A*q1+B*1; % Approx deriv at intermediate value.
17 qstar(:,i+1) = qstar(:,i) + k2*h; % Approx solution at next value of q
18end
19plot(t,yexact,t,qstar(1,:)); % ystar = first row of qstar
20legend('Exact','Approximate');
شکل زیر، جواب دقیق و تقریبی معادله را به ازای $$h$$های مختلف نشان میدهد.
ریاضیات روش رانگ کوتای مرتبه دوم
الگوریتم رانگ کوتای مرتبه دوم به صورت موردی در بالا توضیح داده شد. دیدیم که با استفاده از تقریب مشتق در نقطه میانی بازه بین $$ t_ 0 $$ و $$ t _ 0 + h $$، یعنی در $$ t _ 0 + h / 2 $$، تقریب خوبی برای تابع در $$ t_ 0 + h $$ خواهیم داشت. این تقریب نسبت به حالتی که از مشتق در $$ t _ 0 $$ استفاده کنیم (همان روش اویلر یا رانگ کوتای مرتبه اول) بهتر است. اما آیا راهی برای به دست آوردن روش رانگ کوتای مرتبه دوم از قواعد پایه و اساسی وجود دارد؟ اگر چنین باشد، میتوانیم به تقریبهای بهتری برسیم. در ادامه، این موضوع را بررسی خواهیم کرد.
مرور برخی از ابزارهای ریاضی مفید
در این بخش، دو ابزار مهم ریاضی را معرفی میکنیم که از آنها استفاده خواهیم کرد.
(۱) اگر تابع دومتغیره $$ g ( x , y ) $$ را داشته باشیم، سری تیلور حول $$ (x _0 , y _ 0) $$ به صورت زیر خواهد بود:
$$ \large \eqalign {
g ( { x _ 0 } + \Delta x , { y _ 0 } + \Delta y ) & = g ( { x _ 0 } , { y _ 0 } ) + { \left . { { { \partial g ( x , y ) } \over { \partial x } } } \right | _ { { x _ o } , { y _ 0 } } } \Delta x + { \left . { { { \partial g ( x , y ) } \over { \partial y } } } \right | _ { { x _ o } , { y _ 0 } } } \Delta y + \cdots \cr
& = g + { g _ x } \Delta x + { g _ y } \Delta y + \cdots \cr } $$
که در آن، $$ \Delta x $$ نمو بُعد اول و $$ \Delta y $$ نمو بُعد دوم است. در خط آخر فرمول، از نمادگذاری خلاصه استفاده شده است. که در آن، $$ g _ x $$ نمایانگر مشتق جزئی $$ g $$ نسبت به $$ x $$ و $$ g _ y $$ نشان دهنده مشتق جزئی $$ g $$ نسبت به $$ y $$ است.
(۲) اگر $$ z $$ تابعی از دو متغیر باشد ($$ z = g ( x , y )$$) که در آن، $$ x = r ( t) $$ و $$ y = s ( t) $$، آنگاه با استفاده از قاعده زنجیرهای برای مشتق جزئی داریم:
$$ \large { { d z } \over { d t } } = { { \partial g } \over { \partial r } } { { d r } \over { d t } } + { { \partial g } \over { \partial s } } { { d s } \over { d t } } $$
در حالت خاص، اگر داشته باشیم:
$$ \large y' ( t ) = { { d y ( t ) } \over { d t } } = f ( y ( t ) , t ) = f $$
آنگاه:
$$ \large \begin {align*} y ^ {\prime \prime } ( t ) & = { {d f ( y ( t ) , t ) } \over { d t } } = { { \partial f } \over { \partial t } } { { d t } \over { d t } } + { { \partial f } \over { \partial y } } {{ d y } \over { d t } } \\ & = { { \partial f } \over { \partial t } } + { { \partial f } \over { \partial y } } { { d y } \over { d t } } = { f _ t } + { f _ y } f \end {align*} $$
روش رانگ کوتای مرتبه دوم
برای شروع، باز هم یک معادله دیفرانسیل مرتبه اول را در نظر میگیریم:
$$ \large { { d y ( t ) } \over { d t } } = y' \left ( t \right ) = f ( y ( t ) , t ) $$
اکنون دو تقریب را برای مشتق داریم:
$$ \large \eqalign {
& { k _ 1 } = f (y ^ * ( { t _ 0 } ) , { t _ 0 } ) \cr
& { k _ 2 } = f ( y ^ * ( { t _ 0 } ) + \beta h { k _ 1 } , { t _ 0 } + \alpha h ) \cr } $$
در همه حالات، $$ \alpha $$ و $$ \beta $$ مقادیری کسری بین $$ 0$$ و $$ 1 $$ هستند. این معادلات بیان میکنند که $$ k _ 1 $$ تقریب مشتق بر اساس تخمین مقدار $$ y ( t) $$ در $$ t = t _ 0 $$ (یعنی $$ y ^ * ( t_ 0 ) $$) و در زمان $$ t = t _ 0 $$ است. مقدار $$ k _ 2 $$ مبتنی بر مقدار تخمین زده شده $$ y ^ * ( t _ 0 ) $$ به علاوه کسری از اندازه گام ($$ \beta h $$) ضرب در شیب $$ k _ 1$$ و زمان در $$ t _ 0 + \alpha h $$ (یعنی زمانی بین $$ t _ 0 $$ و $$ t _ 0 + h $$) است. در روشی که در بالا معرفی کردیم، $$ \alpha = \beta = 1 / 2 $$ است، اما انتخابهای دیگری نیز میتوان در نظر گرفت.
برای بهروزرسانی جواب با تقریب بعدی $$ y ( t) $$ در $$ t _ 0 + h $$، از معادله زیر استفاده میکنیم:
$$ \large y ^ * ( { t _ 0 } + h ) = y ^ * ( { t _ 0 } ) + h ( a { k _ 1 } + b { k _ 2 } ) $$
این معادله بیان میکند که مقدار $$ y ^ * ( t _ 0 + h ) $$ از مقدار $$ y ^ * ( t_ 0 ) $$ به علاوه گام زمانی $$h$$ ضرب در یک شیب به دست آمده که مجموع وزنداری از $$ k _ 1 $$ و $$ k_ 2 $$ است.
در روشی که در بالا بیان کردیم، $$ a = 0 $$ و $$ b = 1 $$ است. بنابراین، فقط از تخمین دوم برای شیب استفاده کردیم. لازم به ذکر است که روش اویلر (رانگ کوتای مرتبه اول) حالت خاصی از این روش با $$ a = 1 $$ و $$ b = 0 $$ است و مهم نیست $$ \alpha $$ و $$ \beta $$ چه مقداری داشته باشند، زیرا $$ k _ 2$$ در معادله بهروزرسانی مورد استفاده قرار نمیگیرد.
اکنون هدفمان این است که با استفاده از آن قاعده اول ریاضی که بیان کردیم، تعیین کنیم که چگونه میتوان $$ a $$، $$ b$$، $$ \alpha $$ و $$ \beta $$ را به گونهای یافت که خطا کم شود. از معادله بهروزرسانی بالا شروع میکنیم و میتوانیم عبارات داده شده قبلی را برای $$ k _ 1 $$ و $$ k _ 2 $$ جایگذاری کنیم که منجر به رابطه زیر میشود:
$$ \large { y ^ * } ( { t _ 0 } + h ) = { y ^ * } ( { t _ 0 } ) + h \left ( { a f \left ( { { y ^ * } ( { t _ 0 } ) , { t _ 0 } } \right ) \large + b f \left ( { { y ^ * } ( { t _ 0 } ) + \beta h { k _ 1 } , \,{ t _ 0 } + \alpha h } \right ) } \right )$$
میتوانیم از سری تیلور دوبعدی استفاده کنیم (که نمو بعد اول $$ \beta h k _ 1 $$ و نمو بعد دوم $$ \alpha h $$ است):
$$ \large \eqalign {
f \left ( { { y ^ * } ( { t _ 0 } ) + \beta h { k _ 1 } , \, { t _ 0 } + \alpha h } \right ) & = f + { f _ y } \beta h { k _ 1 } + { f _ t } \alpha h + \cdots \cr
& = f + { f _ y } f \beta h + { f _ t } \alpha h + \cdots \cr } $$
در خط آخر، از این واقعیت استفاده کردهایم که $$ k _ 1 = f $$. جملات مرتبه دوم و بالاتر نوشته نشدهاند. با قرار دادن رابطه اخیر در عبارت قبلی برای $$ y ^ * ( t_ 0 + h )$$ و بازآرایی آن، خواهیم داشت:
$$ \large \eqalign {
{ y ^ * } ( { t _ 0 } + h ) & = { y ^ * } ( { t _ 0 } ) + h \left ( { a f + b \left ( { f + { f _ y } f \beta h + { f _ t } \alpha h + \cdots } \right ) } \right ) \cr
& = { y ^ * } ( { t _ 0 } ) + h a f + h b f + { h ^ 2 } b \beta { f _ y } f + { h ^ 2 } { f _ t } b \alpha + \cdots \cr
& = { y ^ * } ( { t _ 0 } ) + h ( a + b ) f + { h ^ 2 }{ f _ t } b \alpha + { h ^ 2 } b \beta { f _ y } f + \cdots \cr } $$
با استفاده از مشتقات جزئی، میتوان نوشت:
$$ \large \eqalign {
y ( { t _ 0 } + h ) & = y ( { t _ 0 } ) + h { \left . { y' ( t ) } \right | _ { { t _ 0 } } } + { { { h ^ 2 } } \over 2 } { \left . { y ^ {\prime \prime } ( t ) } \right | _ { { t _ 0 } } } + \cdots \cr
& = y ( { t _ 0 } ) + h f + { { { h ^ 2 } } \over 2 } \left ( { { f _ t } + { f _ y } f } \right ) + \cdots \cr
& = y ( { t _ 0 } ) + h f + { { { h ^ 2 } } \over 2 } { f _ t } + { { { h ^ 2 } } \over 2 } { f _ y } f + \cdots \cr } $$
با مقایسه این عبارت با عبارت نهایی تقریبِ
$$ \large { y ^ * } ( { t _ 0 } + h ) = { y ^ * } ( { t _ 0 } ) + h f ( a + b ) + { h ^ 2 } { f _ t } b \alpha + { h ^ 2 } b \beta { f _ y } f + \cdots $$
ثابتهای $$a $$، $$ b$$، $$ \alpha $$ و $$ \beta $$ به صورت زیر هستند:
$$ \large \displaylines {
a + b = 1 \cr
b \alpha = { 1 \over 2 } \cr
b \beta = { 1 \over 2 } \cr } $$
در این دستگاه چهار مجهول و فقط سه معادله وجود دارد؛ بنابراین، بیش از یک جواب خواهیم داشت. واضح است که $$ a = 0$$ و $$ b = 1 $$ و $$ \alpha = \beta = 1/2$$ یکی از جوابها است. از آنجا که همه جملات تقریب برابر با جملات در جواب دقیق هستند، خطای محلی این روش $$ O ( h ^ 3 ) $$ است. خطای سراسری نیز $$ O ( h ^ 2 ) $$ است و به همین دلیل به این روش رانگ کواتای مرتبه دوم میگویند.
فرم دیگر روش رانگ کوتای مرتبه دوم
یک انتخاب دیگر برای ضرایب الگوریتم $$ a = b = 1/ 2 $$ و $$ \alpha = \beta = 1 $$ است. با معادلات $$ k _ 1 $$، $$ k _ 2 $$ و $$ y ^ * ( t _ 0 + h ) $$ شروع میکنیم:
$$ \large \eqalign {
{ k _ 1 } & = f ( { y ^ * } ( { t _ 0 } ) , { t _ 0 } ) \cr
{ k _2 } & = f ( { y ^ * } ( { t _ 0 } ) + \beta h { k _ 1 } , { t _ 0 } + \alpha h ) \cr
{ y ^ * } ( { t _ 0 } + h ) & = { y ^ * } ( { t _ 0 } ) + h ( a { k _ 1 } + b { k _ 2 } ) \cr } $$
با توجه به ضرایبی که گفتیم، داریم:
$$ \large \eqalign {
{ k _ 1 } & = f ( { y ^ * } ( { t _ 0 } ) ,{ t _ 0 } ) \cr
{ k _ 2 } & = f ( { y ^ * } ( { t _ 0 } ) + h { k _ 1 } , { t _ 0 } + \alpha h ) \cr
{ y ^ * } ( { t _ 0 } + h ) & = { y ^ * } ( { t _ 0 } ) + h \left ( { { { { k _ 1 } + { k _ 2 } } \over 2 } } \right ) \cr } $$
بنابرراین، الگوریتم را در این حالت میتوان به صورت زیر نوشت:
تخمین مشتق در $$ t = t _ 0 $$: $$ \large { k _ 1 } = f ( { y ^ * } ( { t _ 0 } ) , \; { t _ 0 } ) $$
تخمین میانی تابع در $$ t = t _ 0 + h $$: $$ \large { y _ 1 } \left ( { { t _ 0 } + h } \right ) = { y ^ * } ( { t _ 0 } ) + { k _ 1 } h $$
تخمین شیب در $$ t = t _ 0 + h $$: $$ \large { k _ 2 } = f \left ( { { y _ 1 } \left ( { { t _ 0 } + h } \right ) , \; { t _ 0 } + h } \right ) $$
تخمین $$ y ( t _ 0 + h )$$ با استفاده از میانگین شیبها: $$ \large { y ^ * } \left ( { { t _ 0 } + h } \right ) = { y ^ * }( { t _ 0 } ) + h \left ( { { { { k _ 1 } + { k _ 2 } } \over 2 } } \right ) $$
برای مثال، باز هم معادله دیفراسیل زیر را در نظر میگیریم:
$$ \large {{dy(t)} \over {dt}} + 2y(t) = 0 $$ یا $$ \large {{dy(t)} \over {dt}} = - 2y(t) $$
که شرایط اولیه آن $$ y ( 0 ) = 3 $$ بوده و میدانیم جواب دقیق آن $$ y ( t) = 3 e ^ {-2t} $$ است. آنچه که میخواهیم انجام دهیم، یافتن مشتق در ابتدای گام زمانی ($$ t = 0 $$ در گام اول) و در پایان گام زمانی $$ t = h = 0.2 $$ است. شیب در نقطه شروع را $$ k _ 1 $$ مینامیم که مقدار عددی آن $$ k _ 1 = - 6 $$ است. شیب در این نقطه در شکل زیر نشان داده شده است. در این شکل، خط (نارنجی) مماس بر منحنی (آبی) در $$ t = 0 $$ مشخص است.
شیب در نقطه انتهایی در شکل زیر نشان داده شده که مقدار عددی آن $$ k _ 2 = -4.0219 $$ است.
حال میانگین دو شیب ($$-5.0110$$) را در نظر گرفته و با استفاده از آن به جلو حرکت میکنیم.
$$ \large {y^*}\left( h \right) = {y^*}(0) + h\left( {{{{k_1} + {k_2}} \over 2}} \right) = 3 + 0.2\left( { - 5.0110} \right) = 1.9978 $$
این جواب به جواب دقیق $$ y ( 0 .2) = 2.0110 $$ بسیاز نزدیک بوده و خطای آن تقریباً $$0.6$$ درصد است. همچنین، این مقدار نزدیک به روش نقطه میانی و بسیار بهتر از روش مرتبه اول است.
الگوریتم روش رانگ کوتای مرتبه دوم (نقطه انتهایی)
معادله دیفرانسیل مرتبه اول زیر داده شده است:
$$ \large {{dy(t)} \over {dt}} = f(y(t),t) $$
برای حرکت از یک نقطه در $$ ( t = t _ 0 , y ^ * ( t_ 0 )) $$ به اندازه یک گام زمانی، مراحل زیر را طی میکنیم:
تقریب مشتق در $$ t = t _ 0 $$: $$ \large { k _ 1 } = f ( { y ^ * } ( { t _ 0 } ) , \; { t _ 0 } ) $$
تقریب مقدار نقطه انتهایی در $$ t = t _ 0 + h $$: $$ \large { y _ 1 } \left ( { { t _ 0 } + h } \right ) = { y ^ * } ( { t_ 0 } ) + { k _ 1 } h $$
تقریب شیب در نقطه انتهایی $$ t = t _ 0 + h $$: $$ \large { k _ 2 } = f \left ( { { y _ 1 } \left ( { {t _ 0 } + h } \right ) , \; { t _ 0 } + h } \right ) $$
تقریب $$ y ( t _ 0 + h ) $$ با استفاده از شیب میانگین: $$ \large { y ^ * } \left ( { { t _ 0 } + h } \right ) = { y ^ * } ( { t _ 0 } ) + h { { \left ( { { k _ 1 } + { k _ 2 } } \right ) } \over 2 } $$
تقریب معادله دیفرانسیل مرتبه سوم با استفاده از متلب
در اینجا، مثال معادله دیفرانسیل مرتبه اول بدون ورودی را برای روش نقطه انتهایی در متلب پیادهسازی میکنیم. در آن مثال از روش نقطه میانی استفاده کردیم. کد متلب مربوط به روش نقطه انتهایی به صورت زیر است:
1%% Example 5
2% Solve y'(t)=-2y(t) with y0=3, endpoint method
3y0 = 3; % Initial Condition
4h=0.2; % Time step
5t = 0:h:2; % t goes from 0 to 2 seconds.
6yexact = 3*exp(-2*t) % Exact solution (in general we won't know this
7ystar = zeros(size(t)); % Preallocate array (good coding practice)
8
9ystar(1) = y0; % Initial condition gives solution at t=0.
10for i=1:(length(t)-1)
11 k1 = -2*ystar(i) % Approx for y gives approx for deriv
12 y1 = ystar(i)+k1*h; % enpoint approximation
13 k2 = -2*y1 % Approx deriv at endpoint.
14 ystar(i+1) = ystar(i) + h*(k1+k2)/2; % Approx solution at next value of y
15end
16plot(t,yexact,t,ystar);
17legend('Exact','Approximate');
با اندکی تغییر در برنامه بالا، اثر اندازه گام $$h$$ بر دقت تقریب قابل بررسی است. شکل زیر نشان میدهد که برای $$h$$های بزرگتر، دقت کاهش پیدا میکند. این نتیجه را برای روش نقطه میانی نیز به دست آوردیم.
اگر این مطلب برایتان مفید بوده است، آموزشهای زیر نیز به شما پیشنهاد میشوند:
- مجموعه آموزشهای محاسبات عددی
- آموزش محاسبات عددی (مرور و حل مساله)
- مجموعه آموزشهای دروس ریاضیات
- آموزش محاسبات عددی با MATLAB
- درون یابی اسپلاین — از صفر تا صد
- درون یابی هرمیت — از صفر تا صد
- درون یابی خطی — به زبان ساده
^^