حل عددی معادلات دیفرانسیل معمولی (IVP)

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

[t,Y] = solver(odefun,tspan,y0)

فرم کلی حل عددی معادلات دیفرانسیل به این صورت است که به جای solver در دستور بالا نوع ode (برای مثال ode45) مورد نظر باید ذکر شود.

منظور از odefun تابعی است که در آن مادلات دیفرانسیل مورئ نظر تعریف شده است و به شکل  زیر است. t متغیر مستقل  اسکالر  و y متغیر وابسته به صورت بردار ستونی  است.

dydt = odefun(t,y)

منظور از tspan بازه ای است که معادله در آن حل می شود و y0 شرایط اولیه است.

مثال)

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

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

q دبی حجمی، h ارتفاع، A سطح مقطع و R مقاومت شیر خروجی مخزن است.

برای حل این معادله ابتدا تابع odefun را تشکیل می دهیم.

فرض کنید A=1, R=1 و تغییرات دبی ورودی به صورت سینوسی باشد.

function dy=testode(t,y)
dy=sin(t)-y;

پس از ذخیره کردن نابع ، در خط فرمان مطلب دستور زیر را می نوسییم :

[t,y]=ode45('testode',[0 ,10],0); 

سری فوریه (2)

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

یکی از راه های بدست آوردن سری فوریه یک تابع استفاده از دستور fit است.این دستور توانایی محاسبه سری فوریه با حداکثر 8 جمله را دارد.

x=[-pi:.1:pi]';

y=sin(x);

f=fit(x,y,'fourier1')

به جای fourier1 می توان  fourier2...fourier8 را قرار داد.

توابع زیر هم نیاز معرفی ندارند/

تابع گاما

Y = gamma(A)

تابع خطا

Y = erf(X)

سری فوریه

فرض کنید تابع f به صورت زیر تعریف شده یاشد.

تابع f یک تابع تناوبی با دوره 2L می باشد و  ثابتهای بسط فوریه  این تابع از روایط زیر بدست می آیند

برای محاسبه این ثوابت از دستور int استفاده می کنیم

R = int(S)
R = int(S,v)
R = int(S,a,b)
R = int(S,v,a,b)

استفاده از این دستور ببسیار ساده است.S عبارتی است که می خواهیم از آن انتگرال بگیریم، v متغیر مستقل است و a,b حدود انتگرال گیری هستند

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

syms w A t n

و بعد از آن محاسبه ثوابت

a0= w/pi*int('A*sin(w*t)','t',0,pi/w)

an=w/pi*int('A*sin(w*t)*cos(n*w*t)','t',0,pi/w)

bn=w/pi*int('A*sin(w*t)*sin(n*w*t)','t',0,pi/w)

قدم بعدی یافتن مقدار ثوابت در n های مختلف است.برای این کار از دستور subs استفاده می کنیم

subs(bn,n,3)

دستور بالا مقدار bn را به ازای n=3 محاسبه می کند

 مقدار an را در n=1  مبهم است  چون صورت و  مخرج an صفر می شود .با استفاده از دستور limit حد این عبارت را در n=1 می توان بدست آورد.

limit(an,n,1,'right')

limit(an,n,1,'left')

تبدیل لاپلاس

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

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

syms t

البته هر متغیر دلخواه دیگری را هم می شود تعریف کرد و تابع برگشتی هم بر حسب s می باشد.

r = laplace(t)

 برای لاپلاس معکوس گرفتن هم از دستور ilaplace استفاده می شود.

ilaplace(r)

حالا فرض کنید بخواهید از تابع 1 لاپلاس بگیرید و بنویسید 

laplace(1)

 نوشتن این دستور به این شکل اشتباه است و موجب خطا می شود، برای این کار باید 1 را به صورت نماد تعریف کنیم

m=sym(1)

ilaplace(m)

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

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

  a=sqrt(sym(2))

f = sym('x^2 + x + 1')

ezplot(f)

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

dsolve('Dy=x')

ans =

x*t+C1

حتما می دانید که عملگر D هم ارز با dy/dx است و هدف ازدستور بالا یافتن جواب معادله dy/dx=x ؟!

متغیر پیش فرض این دستور t است بنابراین x در دستور بالا به صورت یک ثابت عمل می کند.

dsolve('Dy=x','x')

ans =

1/2*x^2+C1

پس هنگامی که متغیر مستقل معادله، t نباشد حتما باید ذکر شود.

در بسیاری از موارد جواب خصوصی معادله مورد نیاز است؛ فرض کنید بخواهیم معادله D2y=x را حل کنیم بطوریکه

y(1)=3

y(2)=5

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

dsolve('D2y=x','y(1)=3,y(2)=5', 'x')

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

 r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')
r = dsolve('eq1','eq2',...,'cond1','cond2',...,'v')

منظور از v متغیر مستقل، cond شرط اولیه و eq معادله دفرانسیل مورد نظر است.

فزض کنید بخواهیم دستگاه معادلات دیفرانسیل زیر را حل کنیم(متغیر مستقل را t فرض کنید).

Dx = y

Dy=-x

به این صورت عمل می کنیم

r=dsolve('Dx = y', 'Dy = -x')

فکر کنم که در مورد گرافیک خیلی حرف زده‌ام و تقربیا همه مطالب اساسی را گفته‌ام و از این به بعد به دیگر قسمتهای مطلب می پردازم .

مطلب یک نرم افزار مهندسی است و اولین و مهمترین چیزی که یک مهندس باید بلد باشد ریاضیات است. اکثر رشته‌های مهندسی در دوره  کارشناسی حدود 19 واحد ریاضی پاس می کنند ولی خیلی ها بلد نیستند در این ضمینه از مطلب استفاده کنند .از این به بعد قرار است یاد بگیرم که ( من هم خیلی چیزها را بلد نیستم) در حد ریاضیات دوره کارشناسی از مطلب استفاده کنیم.

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

هرگونه استفاده تجاری از مــطالب این سایت بصورت كتاب٬ نشریه٬ وب و ... ممنوع می‌باشد

All right reserved©2005 Amin Bashi