|
معادلات دیفرانسیل معمولی – مسایل مقدار اولیه
Ordinary Differential Equation - initial value problems
(IVPs)
با سلام خدمت تمامی دوستان بخصوص انشجویان فیزیک. قصد دارم بتدریج به یکسری از مسایل برنامه نویسی عمومی در زمینه فیزیک که غالبا با آن سرو کار داریم بپردازم. اولین گام ( به شرطی که با مقدمات برنامه نویسی در Matlab خصوصا کار با توابع، ماتریسها آشنا باشید) میتواند حل معادلات دیفرانسیل معمولی باشد.
معادلات دیفرانسیل معمولی به آندسته از معادلات دیفرانسیل اطلاق می شود که تمامی متغیرهای وابسته ( بعنوان مثال x، y و z ) تنها تابعی از یک متغیر مستقل (برای مثال t) می باشند. مثال آشنا در این مورد میتواند معادله نوسانگر هارمونیک باشد که در آن مکان نوسانگر تنها تابعی از زمان است. در مقابل معادله دیفرانسیل معمولی معادلات دیفرانسیل جزیی قرار می گیرند که نمونه بسیار آشنای آن معادله شرودینگر است که در آن تابع موج (متغیر وابسته) تابعی از مختصات مکانی و زمان (یعنی متغیرهای مستقل ) می باشد.
ابتدا به حل مسایل مقدار اولیه (IVP) با معرفی سه روش اویلر (Euler’s method)، رانگ کوتا (Runge–Kutta method) و تابع آماده در متلب بانام ode45 می پردازیم.
1- روش اویلر
این ساده ترین روش است. با استفاده از تقریب مشتق عددی پیشرو الگوریتم تکرار بسادگی بدست می آید :
که در آن h طول گام است. کافست با توجه به معادله دیفرانسیل تابع f(x,t) مشخص گردد و در رابطه تکرار بالا جایگذاری شود. مثال ساده زیر را در نظر بگیرید.
dy/dt+5y=6 , y(t=0)=2
با توجه به معادله بالا خواهیم داشت : f(x,t)=dy/dt=6-5y
بنابراین رابطه تکرار بصورت زیر خواهد بود،
نهایتا برنامه کوتاه زیر را میتوان نوشت :
|
clc,clear all
h=.01;
t0=0;
tf=2;
t=t0:h:tf;
y(1)=2;
N=length(t);
for i=1:N-1
y(i+1)=y(i)+h*(6-5*y(i));
end
plot(t,y)
|
باید توجه داشته باشیم که هر سه روش بکار گرفته در حل مسایل
مقدار اولیه مبتنی بر حل معادله دیفرانسیل مرتبه یک می باشند در حالیکه بسیاری از معادلات دیفرانسیل مهم در فیزیک مرتبه دو هستند. در چنین مواردی با استفاده از تغییر متغیر هر معادله دیفرانسیل مرتبه دو به دو معادله دیفرانسیل مرتبه یک شکافته می شود و هر کدام از این دو معادله بصورت جداگانه با توجه به شرایط اولیه با استفاده از روش اویلر حل میشود. بعنان مثال :
نوسانگر هارمونک ساده
با این توضیحات میتوانیم به حل یک مسئله استاندارد یعنی معادله حرکت نوسانگر هارمونیک بپردازیم. جسمی با جرم معین متصل به فنر تحت تاثیر نیروی بازگردانده در یک بعد حرکت نوسانی انجام می دهد با صرفنظر از نیروهای اتلافی (این مورد در روش رانگ کوتا بررسی می شود) معادله حرکت بشکل زیر است :
با تغییر متغیر v=dx/dt جفت روابط زیر خواهیم داشت،
و بنابراین روابط تکرار مطابق مثال قبلی بدست می آیند،
اگر فرض کنیم که نوسانگر در لحظه t=0 بیشترین فاصله از نقطه تعادل را دارد می توانیم شرایط اولیه زیر را در نظر بگیریم،
برنامه زیر (Oscilator1) معادله نوسانگر هارمونیک ساده را به روش اویلر حل می کند. نمودارها x و v برحسب زمان رسم میشوند و در پایان حرکت واقعی نوسانگر نیز بصورت انیمیشن نمایش داده میشود.
|
function Oscilator1
% simple harmonic Oscilator- Euler's method
% by saeed babanezhad - solid state physics
close all
k=1;
m=1;
w=sqrt(k/m);
h=.1; % time step
t0=0; % initial time
tf=20; % final time
t=t0:h:tf; % time interval
N=length(t);
x=zeros(1,N); % specify size of x matrix
v=zeros(1,N); % specify size of v matrix
x(1)=1; % initial condition for position
v(1)=0; % initial condition for speed
for i=1:N-1
x(i+1)=x(i)+h*v(i);
v(i+1)=v(i)-h*w^2*x(i);
end
figure(1)
subplot(3,1,1)
plot(t,x)
ylabel('x')
subplot(3,1,2)
plot(t,v)
ylabel('y')
% animation
subplot(3,1,3)
% difine spring
spring=line('lineStyle','none','Marker','d',...
'MarkerSize',15,'Color','b','Xdata',[],'ydata',[]);
% difine body
body=line('Marker','s','MarkerSize',30,...
'MarkerFaceColor','r','xdata',[],'ydata',[]);
xlim([min(x) max(x)])
for i=1:N
set(body,'XData',x(i),'YData',0)
set(spring,'XData',linspace(min(x),x(i),20)...
,'YData',zeros(1,20))
drawnow
end
|
در مقابل سادگی روش اویلر خطای آن نسبت به سایر روشها بیشتر است ( از مرتبه 2). در بخش بعدی می خواهم به روش رانگ کوتا ی مرتبه چهار بپردازم (خطا از مرتبه 4) که بسیار دقیقتر از روش اویلر است. |
