حل معادلات دیفرانسیل در نرم افزار متلب با دستور ode45

در این پست قصد داریم ابتدا معادله دیفرانسیل مرتبه اول، سپس معادله دیفرانسیل مرتبه دوم و بالاتر، و در انتها دستگاه معادلات دیفرانسیل را به کمک نرم افزار متلب و کد ode45( دستور ode45) حل کنیم. شکل کلی معادلات دیفرانسیل ODE به صورت زیر است:

${a_{0}}_{(x)}y + {a_{1}}_{(x)}\frac{dy}{dx} + {a_{2}}_{(x)}\frac{d^{2}y}{dx^{2}}+ ...+{a_{n}}_{(x)}\frac{d^{n}y}{dx^{n}} + {b_{0}}_{(x)} = 0$

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

برای شروع ما از حل معادله دیفرانسیل مرتبه اول شروع می کنیم. شکل معادله دیفرانسیل مرتبه اول به صورت $a_{(x)}y + b_{(x)}\frac{dy}{dx} = c_{(x)}$ یا به طور متداول در بین کتاب های درسی به این شکل $y^{'}+ p_{(x)}y =q_{(x)}$ می‌باشد. در پست “معادلات دیفرانسیل و ویژگی های آن” مفاهیم کلی درجه، مرتبه، معادله همگن و ناهمگن و … توضیح داده شده است.

مثال اول: معادله دیفرانسیل مرتبه اول ODE زیر را به همراه شرط اولیه داده شده حل کنید.

$y^{'}=xy$

$x\in[0 ,3] , y_{(x=0)}=1.5$

می خواهیم معادله بالا را به همراه شرایط نوشته شده با کد ode45 حل کنیم. ساختار کلی کد ode45 که ما قصد استفاده از آن را داریم به صورت زیر است:

[t,y] = ode45(odefun,tspan,y0)

همان طور که مشاهده می کنید برای استفاده از این کد، باید معادله دیفرانسیل (odefun) ، بازه‌ی حل معادله (tspan) و شرط اولیه معادله(y0) را به صورت بخصوصی تعریف کنیم. البته شکل کلی تر و توضیحات بیشتر کد ode45 در بانک اطلاعات متلب یعنی اینجا آمده است. در این جا منظور از t متغیر مستقل و منظور از y متغییر وابسته می باشد. در مثال اول و دوم ما، x متغییر مستقل و y متغییر وابسته می باشد.

نحوه تعریف مقدار اولیه y0 در دستور ode45

برای حل مسائل مقدار اولیه با ode45 مقدار تابع در ابتدای بازه‌ی حل را به صورت زیر تعریف می کنیم. اگر مقدار اولیه را نداشته باشیم مسئله از نوع شرط مرزی می باشد و باید با کد bvp4c مسئله را حل کنیم.

y0=1.5

نحوه تعریف بازه حل (tspan) در دستور ode45

از آنجایی که حل ode45 به صورت عددی می باشد، باید بازه‌ای را برای حل تعیین کنیم که جواب ها در آن بازه برای ما به نمایش در بیاید. برای اینکار به صورت زیر بازه‌ی حل را تعریف می کنیم.

1	xspan = [0 3]

در کد ode45 که در بالا آوردیم بازه‌ی حل را به صورت tspan نوشته است زیرا به صورت کلی متغییر مستقل را زمان (time) گرفته است و اصلا اهمیتی ندارد که ما tspan یا xspan یا هر اسم دیگری مثلا baze بنویسم. صرفا جهت زیبایی و شباهت به کد اصلی ما از اسم xspan استفاده کردیم.(چون متغییر مستقل ما x می‌باشد.)

نحوه نوشتن معادله دیفرانسیل(odefun) در دستور ode45

ابتدا باید روی کاغذ معادله را به شکل $y^{'}=F_{(x,y)}$ در بیاوریم. اگر به هیچ وجه این امکان وجود نداشت، یعنی معادله به گونه ای غیر خطی بود که نتوانستیم این کار را بکنیم، باید از روش های دیگر که برای معادلات دیفرانسیل غیر خطی وجود دارد استفاده کنیم.( می توان از روش های خطی سازی معادلات که در CFD وحود دارد استفاده کرد). بعد از تبدیل معادله به فرم گفته شده، آن را به صورت تابع (Function)در نرم افزار متلب می نویسیم.

1
2
3
4
5

function dydx=odefun(x,y)

dydx=y*x

end

عبارت بالا یعنی تابع odefun می باشد و این تابع می تواند دو متغیر x و y را به عنوان ورودی دریافت کند و تابع odefun در ادامه‌ی کد، توسط متغییر dydx تعریف شده است. مثلا اگر به ورودی تابع x=2 و y=3 را بدهیم، خروجی تابع y*x=3*2=6 می‌شود. به طور خلاصه و ساده می‌توان گفت هر جا از عبارت odefun استفاده شود، بیانگر تابع y*x خواهد بود.

در نهایت اگر تمام کد های بالا را کنار هم قرار دهیم، کد نهایی ما برای حل معادله $y^{'}=xy$ به صورت زیر خواهد بود. برای دیدن نتایج در انتهای کد از دستور plot برای رسم نتایج استفاده می کنیم.

1
2
3
4
5
6
7
8
9
10
11
12
13

function example1

xspan = [0 3] ;
y0=1.5 ;
[x,y] = ode45(@odefun,xspan,y0);

plot(x,y)

function dydx=odefun(x,y)
dydx=x*y ;
end

end

فرقی نمی کند که ما تمام کد را به صورت تابع ذخیر کنیم یا غیر تابع. شکل کلی کد به صورت غیر تابع نیز به صورت زیر است:

1
2
3
4
5
6
7
8
9
10
11
12

clear
clc

xspan = [0 3] ;
y0=1.5 ;
[x,y] = ode45(@odefun,xspan,y0);

plot(x,y)

function dydx=odefun(x,y)
dydx=x*y ;
end

حل معادله مرتبه دوم و بالاتر به کمک دستور ODE45

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

$y^{''}+8y^{'}+2y-cos(x)=0$

$x\in[0 ,2] , y_{(x=0)}=0 , y^{'}_{(x=0)}=1$

برای اینکار باید روی کاغذ معادله را به کمک تغییر متغییر به دو معادله ODE تبدیل کنیم. روند تبدیل این صورت می باشد که دو تابع جدید تعریف می کنیم.تایع $y_{1}=y$ و تابع $y_{2}=y^{'}$ .

$y^{'}=y_{2}\rightarrow y^{''}=y_{2}^{'}$

$y=y_{1}\rightarrow y^{'}=y_{1}^{'}=y_{2}$

پس از جاگذاری موارد بالا در معادله اصلی خواهیم داشت:

$\Rightarrow$

$y^{'}_{1}=y_{2}$

$y^{'}_{2}=-8y_{2}-2y_{1}+cos(x)$

بعد از تبدیل معادله دیفرانسیل مرتبه دوم به دو معادله ODE روند حل به صورت قبل می باشد.

نحوه تعریف مقدار اولیه y0 در دستور ode45 برای معادله مرتبه دوم

همانند قسمت قبل می باشد، فقط با این تفاوت که علاوه بر $y_{(x=0)}$ باید مقدار $y^{'}_{(x=0)}$ را هم تعریف کنیم. عدد اول سمت چپ نشانگر $y_{(x=0)}$ و عدد دوم یعنی عدد سمت راست نشانگر $y^{'}_{(x=0)}$ می باشد. ( اگر معادله مرتبه سوم بود، عدد بعدی نشانگر $y^{''}_{(x=0)}$ می شد.)

y0=[0,1]

نحوه تعریف بازه حل (xspan) در دستور ode45 برای معادله مرتبه دوم

دقیقا همانند قسمت قبل می باشد.

1	xspan = [0 2]

نحوه نوشتن معادله دیفرانسیل(odefun) در دستور ode45 برای معادله مرتبه دوم

1
2
3
4
5
6
7
8

function dy=odefun(x,y)

n=2;
dy=zeros(n,1) ;
dy(1)=y(2);
dy(2)=-8*y(2)-2*y(1)+cos(x)

end

عبارت zeros (n,1) به این جهت آورده شده است که باید برای تعریف معادله تابع، بعد از تبدیل معادله مرتبه دوم به دو معادله مرتبه اول ، آن ها را به صورت سطری در زیر هم بنویسیم و تعریف کنیم.به همین دلیل هم هست که خروجی که از تابع ode45 می گیریم دارای دو قسمت است. یعنی تابع y خروجی و نهایی ما شامل دو ستون است که ستون اول آن میزان $y$ و ستون دوم آن میزان $y^{'}$ را نشان می دهد.

در نهایت کد ما به صورت زیر می باشد.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

function example2

xspan = [0 2] ;
y0=[0,1] ;
[x,y] = ode45(@odefun,xspan,y0);
plot(x,y(:,1))

function dy=odefun(x,y)
n=2;
dy=zeros(n,1) ;
dy(1)=y(2);
dy(2)=-8*y(2)-2*y(1)+cos(x);

end

end

برای معادلات مرتبه بالاتر هم به همین صورت هست. به عنوان مثال برای معادله مرتبه ۳ باید با تغییر متغیر $y_{1}=y$ و $y_{2}=y^{'}$ و $y_{3}=y^{''}$ به سه معادله ode برسیم.

حل دستگاه معادلات به کمک دستور ODE45

گاهی اوقات در مسائل، به چند معادله دیفرانسیل وابسته به هم برخورد می کنیم که برای یافتن جواب مسئله، می بایست همزمان یک دستگاه معادلات را حل کنیم. با استفاده از دستور ode45 ما می توانیم دستگاه معادلات دیفرانسیل را به راحتی حل کنیم. روند حل مشابه قسمت قبل یعنی حل معادله دیفرانسیل مرتبه دوم می باشد با این تفاوت که کمی در مرتب سازی معادلات و تغییر متغییر هایی که روی کاغذ انجام می دهیم باید دقت بیشتری بکنیم.

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

$\frac{\text{d}x}{\text{d}t}=−\sigma x + \sigma y$

$\frac{\text{d}y}{\text{d}t}=\rho x - y - xz$

$\frac{\text{d}z}{\text{d}t}=−\beta z + xy$

$\sigma = 10,\beta = \frac{8}{3}, \rho = 28$

$x_{(t=0)} =−۸,y_{(t=0)} = 8, z_{(t=0)} = 27$

برای اینکار باید روی کاغذ معادله را به کمک تغییر متغییر بازنویسی کنیم. روند تبدیل این صورت می باشد. $x_{1}=x$ , $x_{2}=y$ , $x_{3}=z$ و همچنین متغیر مستقل ما زمان هست یعنی t. معادلات بازنویسی شده به صورت زیر می باشند.

$\frac{\text{d}x_{1}}{\text{d}t}=−\sigma x+\sigma y$

$\frac{\text{d}x_{2}}{\text{d}t}=\rho x-y-xz$

$\frac{\text{d}x_{3}}{\text{d}t}=-\beta z+xy$

بعد از بازنویسی دشتگاه معادلات ODE روند حل به صورت قبل می باشد.

نحوه تعریف مقدار اولیه t0 در دستور ode45 برای دستگاه معادلات دیفرانسیل

همانند قسمت قبل می باشد، فقط با این تفاوت که در اینجا x(1)=x , x(2)=y , x(3)=z می باشد.

1	t0=[-8 8 27] ;

نحوه تعریف بازه حل (tspan) در دستور ode45 برای دستگاه معادلات دیفرانسیل

دقیقا همانند قسمت قبل می باشد.

1	tspan = [0 20] ;

نحوه نوشتن معادله دیفرانسیل(odefun) در دستور ode45 برای دستگاه معادلات دیفرانسیل

1
2
3
4
5
6

function dx=odefun(t,x)
n=3;
dx= zeros(n,1);
dx(1)=-10*x(1)+10*x(2);
dx(2)=28*x(1)-x(2)-x(1)*x(3);
dx(3)=-8/3*x(3)+x(2)*x(1);

عبارت zeros (n,1) به این جهت آورده شده است که باید برای تعریف معادله تابع، بعد از تغییر متغیر و بازنویسی دستگاه معادلات ، آن ها را به صورت سطری در زیر هم بنویسیم و تعریف کنیم و در اینجا n تعداد معادلات می باشد(در قسمت قبل n درجه واکنش بود که در نهایت به تعداد درجه واکنش ما به تعداد n معادله می رسیدیم) .به همین دلیل هم هست که خروجی که از تابع ode45 می گیریم دارای سه قسمت است. یعنی تابع x خروجی و نهایی ما شامل سه ستون است که ستون اول آن میزان x(:,1)=x و ستون دوم آن میزان x(:,2)=y و ستون سوم x(:,3)=z را نشان می دهد.

در نهایت کد ما به صورت زیر می باشد.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

function example3

tspan = [0 20] ;
t0=[-8 8 27] ;
[t,x] = ode45(@odefun,tspan,t0);
plot(x(:,1),x(:,3))

function dt=odefun(t,x)

n=3;
dt = zeros(n,1);
dt(1)=-10*x(1)+10*x(2);
dt(2)=28*x(1)-x(2)-x(1)*x(3);
dt(3)=-8/3*x(3)+x(2)*x(1);

end

end

این معادلات معروف به دستگاه معادلات لورنز LORENZ می باشند که شکل های جالبی از رسم آن ها به دست می آید. در زیر نمودار های x-y و x-z و y-z رسم شده اند.

به صورت انیمیشن معادلات lorenz به صورت زیر می باشند.