روش اویلر — از صفر تا صد

در آموزش‌های قبلی از مجموعه مباحث ریاضیات مجله فرادرس، با معادلات دیفرانسیل آشنا شدیم و درباره برخی از روش‌های تحلیلی حل انواع مختلف آن‌ها بحث کردیم. علاوه بر روش‌های تحلیلی، با استفاده از روش‌های عددی نیز می‌توان جواب معادلات دیفرانسیل را محاسبه کرد. روش اویلر، روشی عددی برای حل معادلات دیفرانسیل است. در این آموزش، روش اویلر (Euler Method) یا رونگه-کوتای مرتبه اول (First Order Runge–Kutta) را معرفی می‌کنیم.

روش اویلر

روش اویلر یا رونگه کوتای مرتبه اول، ریاضیات بسیار ساده‌ای دارد و می‌توان آن را به سادگی درک کرد. فرض کنید معادله دیفرانسیل مرتبه اول زیر و شرایط اولیه معلوم $$ y (t _ 0 ) = y _ 0 $$ را داریم:

$$ \large { { d y ( t ) } \over { d t } } = y ^ \prime \left ( t \right ) = f ( y ( t ) , t ) $$

فیلم آموزش معادلات دیفرانسیل به همراه حل سوالات آزمون کارشناسی ارشد در فرادرس

کلیک کنید

تقریب مشتق را به صورت زیر می‌نویسیم:

$$ \large k _ 1 = y ^ \prime \left ( t \right ) = f ( y ^ * ( t ) , t ) $$

بسط تیلور $$ y (t) $$ را حول $$ t _ 0$$ و با فرض گام زمانی $$h$$ می‌نویسیم:

$$ \large y ( { t _ 0 } + h ) = y ( { t _ 0 } ) + y ^ \prime ( {t _ 0 } ) h + y ^ {\prime \prime } ( { t _ 0 } ) { { { h ^ 2 } } \over 2 } + \cdots $$

در عبارت بالا، از جملات بعد از جمله خطی (جملاتی با توان دو و بالاتر) صرفنظر می‌کنیم. بنابراین، عبارت بالا را به صورت زیر می‌نویسیم:

$$ \large y ( { t _ 0 } + h ) \approx y ( { t _ 0 } ) + y ^ \prime ( { t _ 0 } ) h $$

در نتیجه، می‌توانیم جواب تقریبی را در گام زمانی بعدی به صورت زیر به دست آوریم:

$$ \large y ^ * ( { t _ 0 } + h ) = y ^ * ( { t _ 0 } ) + { k _ 1 } h $$

در ادامه، روش اویلر را برای انواع معادلات دیفرانسیل و در قالب مثال بیان می‌کنیم.

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

فرض کنید می‌خواهیم با استفاده از کامپیوتر، جواب معادله دیفرانسیل زیر را به دست آوریم:

$$ \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 ^ { - 2 t } $$ برای $$ t \ge 0 $$ است. این جواب در شکل زیر رسم شده است. با دانستن جواب دقیق، می‌توانیم جواب تقریبی خود را با آن مقایسه کنیم.

راه‌های مختلفی برای به دست آوردن جواب تقریبی وجود دارد. در اینجا از بسط تیلور $$ y ( t) $$ حول $$t=0$$ استفاده می‌کنیم (در حالت کلی، این بسط حول نقطه $$ t = t  _ 0 $$ است). بسط تیلور تابع به صورت زیر است:‌

$$ \large y ( t ) = y ( 0 ) + y ^ \prime ( 0 ) t + y ^ { \prime \prime } ( 0 ) { { { t ^ 2 } } \over 2 } + \cdots $$

اکنون جواب را به گام زمانی کوچک $$h$$ بعد از $$t =0 $$ محدود کرده و بسط تیلور را تا مشتق اول در نظر می‌گیریم. بنابراین، داریم:

$$ \large \eqalign {
& y ( h ) = y ( 0 ) + y^ \prime ( 0 ) h + y ^ {\prime \prime } ( 0 ){ { { h ^ 2 } } \over 2 } + \cdots \cr
& y ( h ) \approx y ( 0 ) + y ^ \prime ( 0 ) h \cr } $$

در ادامه، تقریب جواب را $$ y ^ * (h) $$ و مشتق آن را $$ y ^ \prime ( 0 ) =  k _ 1 $$ می‌نامیم.

$$ \large y ^ * ( h ) = y ( 0 ) + { k _ 1 } h $$

شکل زیر، نمودار جواب را برای $$ h = 0.2 $$ نشان می‌دهد.

نمودار سمت راست، ناحیه سایه زده شده نمودار سمت چپ را با بزرگنمایی نشان می‌دهد. تقریب تیلور مرتبه اول، یک خط راست است که از مقدار اولیه شروع شده و شیب آن $$-6$$ است (یعنی $$ k _ 1 = y ^ \prime ( 0 ) = -2 y (0) = - 6 $$).

تقریب در $$ t= h = 0.2 $$ برابر با حاصل‌ جمع مقدار اولیه و ضرب شیب در گام زمانی $$h$$ است:

$$ \large y ^ * ( h ) = y ^ * ( 0 .2 ) = y ( 0 ) + { k _ 1} y ^ \prime (0) = 1.8$$

فیلم آموزش حل معادلات دیفرانسیل با متمتیکا Mathematica در فرادرس

کلیک کنید

واضح است که هرچه $$h$$ کوچکتر باشد، تقریب بهتر و دقیق‌تر خواهد بود.

برای پیدا کردن تقریب گام بعدی $$ y ^ * ( 2 h ) $$، به سادگی فرایند را تکرار کرده و از مقدار قبلی $$ y ^ * ( h ) $$ برای تخمین مشتق در زمان $$h$$ استفاده می‌کنیم (مقدار دقیق $$ y ( h ) $$ را نمی‌دانیم؛ بنابراین، فقط می‌توانیم مشتق ($$ k _ 1 $$) را محاسبه کنیم). ترتیب انجام مراحل به صورت زیر است:

عبارت دقیق مشتق:     $$ \large y ^ \prime ( t ) = - 2 y ( t ) $$

تقریب مشتق:    $$ \large { k _ 1 } = - 2 y ^ * ( h ) $$

سری تیلور حول $$t = h $$:     $$ \large y ( 2 h ) = y ( h ) + y ^ \prime ( h ) h + y ^ { \prime \prime }( h ) { { { h ^ 2 } } \over 2 } + \cdots $$

بسط تیلور بریده:     $$ \large y ( 2 h ) \approx y ( h ) + y ^ \prime ( 0 ) h $$

جواب تقریبی:    $$ \large y ^ * ( 2 h ) = y ^ * ( h ) + { k _ 1 } h $$

حرکت به سمت جلو از زمان $$  t _0 $$ تا $$ t_0 + h $$ در حالت کلی به صورت زیر خواهد بود:

عبارت دقیق مشتق در $$ t = t_0$$:     $$ \large y ^ \prime ( t_0 ) = - 2 y ( t_0 ) $$

تقریب مشتق با استفاده از تقریب قبلی $$ y (t ) $$:    $$ \large { k _ 1 } = - 2 y ^ * ( t_0 ) $$

سری تیلور حول $$t = t _ 0 $$:     $$ \large y ( t_0 + h ) = y ( t _ 0 ) + y ^ \prime ( t _ 0 ) h + y ^ { \prime \prime }( t _ 0 ) { { { h ^ 2 } } \over 2 } + \cdots $$

بسط تیلور بریده:     $$ \large y ( t _ 0 + h ) \approx y ( t_ 0) + y ^ \prime ( t_ 0 ) h $$

جواب تقریبی:    $$ \large y ^ * ( t _ 0 +h ) = y ^ * ( t _ 0 ) + { k _ 1 } h $$

پیاده‌سازی در متلب

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

1%% Example 1
2% Solve y'(t)=-2y(t) with y0=3
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);  % Previous approx for y gives approx for derivative
12    ystar(i+1) = ystar(i) + k1*h; % Approximate solution for next value of y
13end
14plot(t,yexact,t,ystar);
15legend('Exact','Approximate');

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

خلاصه الگوریتم

برای به دست آوردن جواب تقریبی معادله دیفرانسیلِ

$$ \large { { d y ( t ) } \over { d t } } = f ( y ( t ) , t ) $$

به روش رونگه-کوتای مرتبه اول، از نقطه $$t _0 $$ شوع کرده و با گام زمانی $$ h $$، مراحل زیر را به ترتیب طی می‌کنیم:

تقریب مشتق:    $$ \large { k _ 1 } = f ( { y ^ * } { \rm { ( } } { { \rm{ t }} _ 0 } ) , { t _0 } ) $$

تقریب جواب در گام زمانی بعدی:    $$ \large { y ^ * } ( { t _ 0 } + h ) = { y ^ * } ( { t _ 0 } ) + { k _ 1 } h $$

داشتن مقدار اولیه برای شروع الگوریتم، ضروری است.

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

اگر معادله مرتبه اول، ورودی داشته باشد، یا اصطلاحاً ناهمگن باشد، مشکلی خاصی برای جواب تقریبی وجود ندارد، اما یافتن جواب دقیق یا واقعی سخت‌تر خواهد شد. مثال زیر را در نظر بگیرید:

$$ \large y ^ \prime ( t ) + 2 y ( t ) = \cos ( 4 t ) $$

فیلم آموزش محاسبات عددی با متلب MATLAB در فرادرس

کلیک کنید

یا

$$ \large y ^ \prime ( t ) = - 2 y ( t ) + \cos ( 4 t ) $$

برای به دست آوردن جواب تقریبی، مراحل زیر را طی می‌کنیم:

عبارت دقیق مشتق در $$ t = t _0 $$:    $$ \large y ^ \prime ( { t _ 0 } ) = - 2 y ( { t _ 0 } ) + \cos ( 4 { t _ 0 } ) $$

تقریب مشتق:    $$ \large { k _ 1 } = - 2 y ^ * ( { t _ 0 } ) + \cos ( 4 { t _ 0 } ) $$

بسط تیلور حول $$ t = t _0 $$:    $$ \large y ( { t _ 0 } + h ) = y ( { t _ 0 } ) + y ^ \prime ( { t _ 0 } ) h + y ^ { \prime \prime } ( { t _ 0 } ) { { { h ^ 2 } } \over 2 } + \cdots $$

بسط تیلور بریده:    $$ \large y ( { t _ 0 } + h ) \approx y ( { t _ 0 } ) + y ^ \prime ( { t _ 0 } ) h $$

جواب تقریبی:    $$ \large y ^ * ( { t _ 0 } + h ) = y ^ * ( { t _ 0 } ) + { k _ 1 } h $$

جواب دقیق این مثال $$ y ( t ) = 2.9 { e ^ { - 2 { \kern 1 pt } t } } + 0.1 \cos ( 4 t ) + 0.2 \sin ( 4 t ) $$ است.

پیاده‌سازی در متلب

برای پیاده‌سازی محاسبات بالا، می‌توانیم از نرم‌افزار متلب استفاده کنیم. دقت کنید که در محاسبات، مقدار $$k_1$$ تغییر می‌کند.

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.
6% Exact solution, hard to find in this case (in general we won't have it)
7yexact = 0.1*cos(4*t)+0.2*sin(4*t)+2.9*exp(-2*t);
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 = -2*ystar(i)+cos(4*t(i));  % Previous approx for y gives approx for deriv
13    ystar(i+1) = ystar(i) + k1*h; % Approximate solution for next value of y
14end
15plot(t,yexact,t,ystar);
16legend('Exact','Approximate');

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

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

حل معادلات دیفرانسیل غیرخطی به صورت تحلیلی، کار دشواری است. اما در تقریب زدن آن، مشکلی وجود ندارد. معادله دیفرانسیل زیر را در نظر بگیرید.

$$ \large { { d y ( t ) } \over { d t } } - y ( t ) ( 1 - 2 t ) = 0 , \quad \quad \quad \quad y ( 0 ) = 1 $$

جواب دقیق این معادله به صورت زیر است:

$$ \large y ( t ) = { e ^ { t - { t ^ 2 } } } $$

برای به دست آوردن جواب تقریبی، عبارت $$ y ^ \prime  (t) $$ را نوشته و از آن برای یافتن $$ k _ 1 $$ استفاده می‌کنیم. محاسبات را به صورت زیر انجام می‌دهیم:

عبارت دقیق مشتق در $$ t = t _ 0 $$:    $$ \large y ^ \prime ( { t _ 0 } ) = y ( { t _ 0 } ) ( 1 - 2 { t _ 0 } ) $$

تقریب مشتق:    $$ \large { k _ 1 } = y ^ * ( { t _ 0 } ) ( 1 - 2 { t _ 0 } ) $$

جواب تقریبی:    $$ \large y ^ * ( { t _ 0 } + h ) = y ^ * ( { t _ 0 } ) + { k _ 1 } 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  ystar(i+1) = ystar(i) + k1*h; % Approximate solution at next value of y
14end
15plot(t,yexact,t,ystar);
16legend('Exact','Approximate');

نمودار جواب دقیق و جواب‌های تقریبی به ازای مقادیر مختلف $$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  \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 } $$

که در آن، $$\gamma (t) $$ تابع پله واحد است.

جواب دقیق این معادله برای $$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 } } $$ است.

با تعریف متغیرهای جدید $$q _ 1$$، $$ q _ 2$$ و $$q _ 3$$ معادله دیفرانسیل فوق را با سه معادله مرتبه اول توصیف می‌کنیم:

$$ \large \eqalign {
{ q _ 1 } ( t ) & = y ( t ) \cr
{ q _ 2 } ( t ) & = y ^ \prime ( t ) \cr
{ q _ 3 } ( t ) & = y ^ { \prime \prime } ( t ) \cr
\cr
{ q _ 1 } ^ \prime ( t ) & = y ^ \prime ( t ) & = { q _ 2 } (t ) \cr
{ q _ 2 } ^ \prime ( t ) & = y ^ { \prime \prime } ( t ) & = { q _ 3 } ( t ) \cr
{ q _ 3 } ^ \prime ( t ) & = y ^ {\prime \prime\prime } ( t ) & = \gamma ( t ) - 4 y ^ { \prime \prime } ( t ) - 6 y ^ \prime ( t ) - 4 y ( t ) \cr
& & = \gamma ( t ) - 4 { q _ 3 } ( t ) - 6 { q _ 2 } ( t ) - 4 { q _ 1 } ( t ) \cr } $$

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

$$ \large \eqalign {
{ \bf { q } } ^ \prime ( t ) & = \left [ { \matrix {
{ { q _ 1 } ^ \prime ( t ) } \cr
{ { q _ 2 } ^ \prime ( t ) } \cr
{ { q _ 3 } ^ \prime ( 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 } } ^ \prime ( 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 } $$

اکنون، مسئله مشابه موارد قبلی است؛ با این تفاوت که متغیرها و مقادیر $$k_1$$، در این حالت بردار هستند:

عبارت دقیق مشتق در $$ t = t_0 $$:    $$ \large { \bf { q } } ^ \prime ( { t _ 0 } ) = { \bf { A q } } ( { t _ 0 } ) + { \bf { B \gamma } } ( { t _ 0 } ) $$

تقریب مشتق:    $$ \large { { \bf { k } } _ 1 } = { \bf { A } } { { \bf { q } } ^ * } ( { t _ 0 } ) + { \bf { B \gamma } } ( { t _ 0 } ) $$

بسط تیلور حول $$ t = t _0$$:    $$ \large { \bf { q } } ( { t _ 0 } + h ) = { \bf { q } } ( { t _ 0 } ) + { \bf { q } } ^ \prime ( { t _ 0 } ) h + { \bf { q } } ^ {\prime \prime } ( { t _ 0 } ) { { { h ^ 2 } } \over 2 } + \cdots $$

بسط تیلور بریده:    $$ \large { \bf { q } } ( { t _ 0} + h ) \approx { \bf { q } } ( { t _ 0 } ) + { \bf { q } } ^ \prime ( { t _ 0 } ) h $$

جواب تقریبی:    $$ \large { \bf { q } } ^ * ( { t _ 0 } + h ) = { \bf { q } } ^ * ( { t _ 0 } ) + { { \bf { k } } _ 1 } h $$

پیاده‌سازی در متلب

مانند مثال‌های قبلی، محاسبات اخیر را می‌توانیم در متلب پیاده‌سازی کنیم. برنامه مربوط به این مثال به صورت زیر است:

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  qstar(:,i+1) = qstar(:,i) + k1*h; % Approximate solution at next value of q
16end
17plot(t,yexact,t,qstar(1,:));        % ystar = first row of qstar
18legend('Exact','Approximate');
19hold off

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