روش رانگ کوتای مرتبه دوم – از صفر تا صد

تعریف مسئله

$$\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$$). در ادامه، یک روش شهودی برای انجام این کار ارائه کرده و مثال‌هایی را نیز بررسی می‌کنیم. این روش همان روش رانگ کوتای مرتبه دوم است.

روش رانگ کوتای مرتبه دوم

در این بخش، روش رانگ کوتای مرتبه دوم را برای دسته‌هایی از معادلات دیفرانسیل به صورت شهودی معرفی خواهیم کرد.

معادله دیفرانسیل مرتبه اول بدون ورودی

$$\large { { d y ( t ) } \over { d t } } + 2 y( t ) = 0$$   یا  $$\large { { d y ( t ) } \over { d t } } = - 2 y ( t )$$

در شکل زیر، شیب در $$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 )$$

سری تیلور بریده:    $$\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$$

تقریب معادله دیفرانسیل مرتبه اول بدون ورودی با استفاده از متلب

هر خط از کد بالا به خوبی مراحل الگوریتم را نشان می‌دهد. با تغییر اندکی در برنامه، می‌توان اثر اندازه $$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$$ را با استفاده از مقدار مناسب متغیر زمان تغییر دهیم. برنامه زیر در متلب مربوط به مثال اخیر است.

با کمی تغییر در برنامه بالا، می‌توانیم تقریب را برای $$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$$

تقریب معادله دیفرانسیل غیرخطی مرتبه اول با ورودی با استفاده از متلب

مانند مثال‌های قبل، به سادگی می‌توانیم برنامه تقریب معادله دیفرانسیل غیرخطی مرتبه اول با ورودی را در متلب بنویسیم. این برنامه به صورت زیر است:

شکل زیر، جواب دقیق و تقریبی معادله را به ازای $$h$$های مختلف نشان می‌دهد.

همان‌طور که مشاهده می‌کنیم، با کاهش $$h$$، دقت تقریب زیاد می‌شود.

معادله دیفرانسیل خطی مرتبه بالاتر

معادله دیفرانسیل مرتبه بالاتر زیر را در نظر بگیرید:

$$ \large \displaylines {<br /> { { { 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<br /> { \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 } $$