تبلیغات
آموزش مطلب - pdepe

pdepe

نویسنده :امین باشی
تاریخ:یکشنبه 5 شهریور 1385-12:08 ب.ظ

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

فرض کنید u تابعی باشد که در معادله زیرصدق کند.

 

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

این معادله می تواند معادله حاکم بر انتقال حرارت در یک تیغه , استوانه و یا کره باشد؛ در این صورت  f عبارت مربوط به شار انتقال حرارت و s مربوط به تولید یا مصرف انرژی می باشد.

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

  sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)

m : مشخص کننده هندسه مساله است. 0 برای تیغه، 1 برای سیلندر و 2 برای کره

pdefun : تابعی که معادله را تعریف می کند

[c,f,s] = pdefun(x,t,u,dudx)

c, f, s همان پارامترهای معادله دیفرانسیل پاره ای هستند

icfun : تابعی که شرایط اولیه را تعریف می کند

u = icfun(x)

bcfun : تابعی که شرایط مرزی را بیان می کند

[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

اندیس l مربوط x0 , و اندیس r مربوط به xn (نقاط ابتدایی و انتهایی بردار xmesh )

xmesh : برداری شامل  نقاط x1 تا xn 

tspan : بردار زمان متناظر با بردار xmesh 

مثال

یک لوله استوانه ای را در نظر بگیرید که از وسط با یک غشا به دو نیم تقسیم شده است.در یک طرف این لوله گاز A با فشار 10 بار و در طرف دیگر گاز B وجود دارد.اگر t=0 غشا پاره شود گاز A در B نفوذ می کند تغییرات فشار جزیی گاز A را در طول لوله حساب کنید.

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

D ضریب نفوذ گاز A در B است و معمولا از مرتبه 1e-5 است.

شرایط اولیه

و شرایط مرزی

با مقایسه معادله3 با معادله 1 می بینیم که

m = 0

c = 1

s = 0

f = 1e-5 * DuDx

پس تابع pdefun به این صورت تعریف می شود.

function [c,f,s] = pdefun0(x,t,u,DuDx)
c = 1;
f = 1e-5*DuDx;
s = 0;

تابع pdeic

function u0 = pdeic0(x)
if ((x >= 0) & (x <= .5))
     u0=10;
elseif ((x >= 0.5) & (x <= 1))
    u0=0;
end

 و تابع pdebc

function [pl,ql,pr,qr] = pdebc0(xl,ul,xr,ur,t)
pl = 0;
ql = 100000;
pr = 0;
qr = 100000;

حالا با استفاده از pdepe (برنامه زیر) می توانیم معادله را حل کنیم.

function pdex
x=linspace(0,1,20);
t=linspace(0,30000,20);
m=0;
sol = pdepe(m,@pdefun0,@pdeic0,@pdebc0,x,t);
u=sol(:,:,1);
uu=u;
surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')
%--------------------------------------------
function [c,f,s] = pdefun0(x,t,u,DuDx)
c = 1;
f = 1e-5*DuDx;
s = 0;
%---------------------------------------------
function u0 = pdeic0(x)
if ((x >= 0) & (x <= .5))
u0=10;
elseif ((x >= 0.5) & (x <= 1))
u0=0;
end
%---------------------------------------------
function [pl,ql,pr,qr] = pdebc0(xl,ul,xr,ur,t)
pl = 0;
ql = 100000;
pr = 0;
qr = 100000;
 








All right reserved©2005 Amin Bashi