حل معادلات دیفرانسیل در MATLAB

معادلات دیفرانسیل معمولی (ODE) به آن دسته از معادلات دیفرانسیل گفته می‌شود که توابع مجهول در آن‌ها بر حسب تنها یک متغیر مستقل باشد. معادله‌ای که بیشتر از یک متغیر مستقل داشته باشد را معادله دیفرانسیل نسبی یا معادله دیفرانسیل با مشتقات جزئی (PDE) می‌نامند.

دستور dsolve برای حل معادله دیفرانسیلی در متلب، به کار می رود. اگر y تابعی از متغیر x باشد و معادله دیفرانسیلی شامل مشتق مرتبه اول و یا مرتبه های بالاتر از y باشد، می توان از دستور dsolve جهت حل معادلات دیفرانسیلی استفاده نمود. در این دستور از حرف کلیدی D به عنوان نماد مشتق استفاده می شود.

مثال: میخواهیم معادله دیفرانسیلی xy′+۱=y را حل کنیم:

1	dsolve('x*Dy+1=y','x')

 نتیجه:

1 2	ans = C2*x + 1

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

همان طور که گفته شد، مشتق مرتبه اول را با نماد Dy نشان می دهیم.در دستور dsolve اگر مشتق مرتبه دوم و یا بالاتر باشد، آنگاه باید ابتدا نماد D را نوشته، سپس عدد مربوط به مرتبه مشتق را بنویسیم و در آخر نیز نماد y نوشته شود. مثلا برای تعریف مشتق دوم تابع y باید D2y و برای تعریف ”’y باید نماد D3y را به کار ببریم.

جهت حل معادلات دیفرانسیلی بصورت پارامتریک می توان پارامترها به سادگی در معدله قرار داد. در این روش نیازی به دستور syms نمی باشد. بطور مثال جهت حل معادله ax”+bx’+cx=d می نویسیم:

1	w=dsolve('a*D2x+b*Dx+c*x=d')

که D2 به معنی مشتق دوم و D به معنی مشتق اول می باشد. جواب معادله فوق در متغیر w ذخیره می شود.

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

1 2	x''+2x'+x=1 w=dsolve('D2x+2*Dx+x=1')

 جواب معادله:

1	w=C20*exp(-t) + C21*t*exp(-t) + 1

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

برای حل معادله y’= y +sin (t) , y(0) =1 به شکل زیر عمل می کنیم:

1 2	w=dsolve('Dy = y + sin(t)' , ' y (0)=1') w=(3*exp(t))/2 - cos(t)/2 - sin(t)/2

در صورتی که بخواهیم جواب بدست آمده را رسم کنیم باید از دستور ezplot استفاده کنیم

1

ezplot(w)

 نتیجه:

حل معادلات دیفرانسیلی سخت و غیرسخت

بسته به نوع معادلات که اصطلاحا به آنها سخت (stiff) و غیرسخت (stiffness) گفته می شود، روش حل آنها در MATLAB کمی متفاوت خواهد بود. اصطلاح سخت (stiff) برای آن دسته از معادلاتی بکار می رود که برای مثال در مقابل متغیر مستقلی همچون t چند متغیر وابسته مانند x و y و… وجود دارد، بگونه ای که اندازه مشتقات متغیرهای وابسته نسبت به متغیر مستقل بطور قابل ملاحظه ای متفاوت است. در غیر اینصورت معادله غیرسخت نامیده می شود. همچنین معادلات سخت شامل آن دسته از معادلات دیفرانسیلی می شوند که حل آنها با روش محاسبات عددی پایدار و همگرا نبوده و تنها راه حل آنها، بسیار کوچک کردن گام (Step) در روش عددی می باشد. اگر در معادله دیفرانسیلی، متغری وجود دارد که باعث تغییرات بسیار زیاد در جواب مساله می شود، این دسته را جز معادلات غیرسخت طبقه بندی می کنند.

کاربرد	دقت	نوع مساله	دستور
اکثر موارد (سعی شود جهت حل معادله ابتدا از این دستور استفاده شود)	متوسط	غیرسخت	ode45
حل مسائل دارای خطای خام (crude error)، حل مسائل تقریبا سخت	پایین	غیرسخت	ode23
حل مساله دارای خطای دقیق، معادلات مربوط به محاسبات عددی زمان بر	پایین تا بالا	غیرسخت	ode113
هنگامی که حل معادله با دستور ode45 بسیار کند پیش رود.	پایین تا متوسط	سخت	ode15s
مسائل دارای خطای خام (crude error) ، هنگامی که ماتریس جرم (mass matrix) ثابت باشد	پایین	سخت	ode23s
حل معادلات بدون میرایی عددی (numerical damping)	پایین	سخت	ode23t
حل مساله با خطای خام و معادله سخت	پایین	سخت	ode23tb

مثال) حل معادله زیر با دستور ode45

گام اول: بازنویسی مسئله بشکل مجموعه ای از معادلات مرتبه یک. با تغییر متغیر y۱‘=y۲ دو معادله دیفرانسیل مرتبه یک حاصل میشود:

گام دوم: استفاده از ODE های مرتبه یک، پس از بازنویسی معادلات مطابق گام اول، باید ان معادلات را برای بکارگیری ode45 کد نماییم:

توجه نمایید که تابع باید دارای دو آرگومان باشد؛ یکی مربوط به متغیر مستقل (t) و دیگری مربوط به متغیرهای وابسته (y). حتی اگر هیچکدام در تابع های مرتبه یک ظاهر نشوند، باید آن را لحاظ نمود:

1 2 3 4 5

function dydt=vdp1(t,y) epsilon=5; w=2.466; f=5; dydt=[y(2) ; epsilon*(1-y(1)^2)*y(2)-y(1)+f*cos(w*t)];

 گام سوم- استفاده از یک حل کننده (solver) برای مسئله، همانگونه که گفته شد، قصد داریم از دستور ode45 استفاده نماییم. تابع ode45 سه ورودی(آرگومان) دارد. اولین آرگومان تابعی است که آن را در گام دوم ایجاد کرده ایم، دومین آرگومان بازه مربوط به متغیر مستقل است که در مورد این مساله بازه زمانی t است و نهایتا آرگومان سوم بردار شرایط اولیه است. در این مثال بازه زمانی [۰ ۱۰۰] و مقادیر اولیه y(1)=1.2 و y(2)=0 را مطابق زیر در ode45 قرار می دهیم:

1	[t,y]=ode45(@vdp1,[0 100],[1.2 0]);

در این مثال برای پذرفتن vdp1 بعنوان یک تابع از نماد @ استفاده شده است. خروجی ode45 برداری از زمان t و ماتریس پاسخ y می باشد.

گام چهارم- مشاهده خروجی، اینکار را براحتی با استفاده از plot می توان انجام داد.

1 2 3

plot(y(:,1),y(:,2)) xlabel('y1') ylabel('y2')

 خروجی:

با این توضیحات برنامه vdp1 با فرم نهایی و ساده زیر، نوسانگر واندپل را در فضای فازی شبیه سازی میکند.

1 2 3 4 5 6 7 8 9 10

function vonerpol [t,y]=ode45(@vdp1,[0 100],[1.2 0]); plot(y(:,1),y(:,2)) xlabel('y1') ylabel('y2') function dydt=vdp1(t,y) epsilon=5; w=2.466; f=5; dydt=[y(2) ; epsilon*(1-y(1)^2)*y(2)-y(1)+f*cos(w*t)];

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

که تابع f بصورت f(t,y)=dy/dt می باشد.

خطای این روش از مرتبه ۴ می باشد.

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

روش اویلر

روش هون

روش تیلور

روش رانگ-کوتا

روش آدامز-بشفورت-مولتون

روش میلن-سیمپسون

روش هامینگ

روش رانگ-کوتا-فلبرگ مرتبه ۵

روش رحمانزاده-کای-وایت

تمام روش‌ها بیان شده راه حل دقیق را ایجاد نمی‌کند و فقط یک تقریب به‌دست می‌آید. چون این روش‌ها دارای محاسبات زیادی هسند، تنها جواب‌هایی در فواصل زمانی مجزا می‌دهند. مشخصا جواب‌ها در زمان ابتدایی شرایط وفاصله زمان‌های مشخص، h، بدست می‌آید. (t۰, t۱ = t۰ + h, t۲ = t۰ + ۲h)