همه ما می توانیم  معادله a*x^2+b*x+c=0 به صورت تحلیلی حل کنیم ولی فکر نمی کنم حل تحلیلی این معادله را بلد باشیم a*x^3+b*x^2+c*x+d=0 البته یک راه تقریبی برای حل آن هست ولی کاربرد زیادی ندارد.با استفاده از مطلب پیدا کردن حل تحلیلی این معادله بسیار ساده است

solve('a*x^3+b*x^2+c*x+d')

دستور solve برای یافتن حل تحلیلی معادات جبری بکار می رود

 g = solve(eq)
g = solve(eq,var)
g = solve(eq1,eq2,...,eqn)
g = solve(eq1,eq2,...,eqn,var1,var2,...,varn)

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

 

g= solve('a+b=1','a')

class(g)

دستور بالا مقدار a را طوری پیدا می کند که عبارت a+b=1 برقرار باشد و دستور زیر این معادله برای b حل می کند

g= solve('a+b=1','b')

وقتی بخواهیم چند معادله را همزمان حل کنیم خروجی  دستور یک structure از نوع سمبولیک است

 g= solve('a^2+b=2','a-b')

class(g)

g.a

g.b

[ شنبه 25 شهریور 1385 - 04:09 ق.ظ ]
[ویرایش شده در : - - -]

[ پیام ()|| امین باشی ] [solve , ] [+]

معمولا دستگاه معادلات  خطی  را به صورت زیر نمایش می دهند

A X = B

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

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

syms a1 a2 a3 a4

syms x1 x2

A=[a1 a2;a3 a4]

X=[x1;x2]

B=A*X

X=A^-1*B

تجزیه LU

A = pascal(3);

[L,U]=lu(A)

روش حذفی گوس

بهترین روش حل دستگاه معادلت خطی (سریع تر و دقیق تر) استفاده از روش حذفی گوس است

اگر دستگاه را به صورت زیر در نظر بگیریم

 A X = B

آنگاه

X = A\B

و اگر دستگاه به این صورت باشد

 XA = B

خواهیم داشت

X = B/A

استفاده از شکل دوم دستگاه چندان مرسوم نیست و کمتر مورد استفاده قرار می گیرد

A=[a1 a2;a3 a4]

X=[x1;x2]

B=A*X

X=A\B

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

فرض کنید که y , t  به صورت زیر موجود باشند

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';

و بخواهیم آنها با صورت

y(t) = c1 + c2 * exp(-t)

برازش کنیم . در این حالت دو مجهول و 6 معادله داریم یعنی دستگاه Overdetermined (معادل فارسی اش را بلد نیستم)است.

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

f=fit(t,y,'c1+c2*exp(-x)')

اما روش دیگری هم برای حل معادله وجود دارد 

ماتریس E که یک متاریس 2*6 است را طوری انتخاب می کنیم که در رابطه زیر صدق کند

EC = y

C ماتریس  مجهولات است

E = [ones(size(t)) exp(-t)]

C = E\y

[ سه شنبه 21 شهریور 1385 - 04:09 ق.ظ ]
[ویرایش شده در : - - -]

[ پیام ()|| امین باشی ] [حل معادلات خطی , ] [+]

حل معادلات BVP

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

sol = bvp4c(odefun,bcfun,solinit)

فرض کنید y در معادله زیر صدق کند

D2y = sin(x) + 2x

و داشته باشیم

y(-pi) = 10.8696

y(pi) = 10.8696

با توجه به شرایط داده شده مساله از نوع BVP  (فرق BVP با IVP رو می دانید؟) و از مرتبه دوم است.هر وقت مرتبه معادله دفرانسیل از یک بیشتر باشد، باید آن را به معادلات مرتبه اول تجزیه کنیم.

اگر

y(1) = y

y(2) = Dy

خواهیم داشت

dydx = [ y(2);-sin(x)+2*x];

حالا شرایط مرزی را می نویسیم :

اگر res مشخص کننده شرایط مرزی باشد، شرایط مرزی را طوری می نویسیم که مقدار res در آن نقاط صفر شود یعنی

res = [ ya(1)-10.8969;yb(1)-10.8969];

منظور از a,b نقاطی است که مقدارمرزی در آن تعریف شده است.

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

function bvptest
solinit = bvpinit(linspace(-pi,pi,10),[-1 0]);
sol = bvp4c(@twoode,@twobc,solinit);
x = linspace(-pi,pi);
y = deval(sol,x);
plot(x,y(1,:));
%----------------------------------------------
function dydx = twoode(x,y)
dydx = [ y(2);-sin(x)+2*x];
%----------------------------------------------
function res = twobc(ya,yb)
res = [ ya(1)-10.8969;yb(1)-10.8969];

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

solinit = bvpinit(x,yinit)

بردار x بازه ای است که معادله در آن حل می شود. و y حدس اولیه است.

در این مثال، معادله را در بازه [pi,pi-] حل کرده و این بازه را به 10 قسمت تقسیم کرده ایم.

کار دستور deval محاسبه جواب معادله در یک بازه دیگر است، این بازه باید زیر مجموعه بازه قبلی باشد.

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

[ شنبه 11 شهریور 1385 - 05:09 ق.ظ ]
[ویرایش شده در : - - -]

[ پیام ()|| امین باشی ] [BVP , ] [+]

تا حالا شده کخ بخوهایم محتویات یک دایرکتوری را چاپ کنید؟واسه من که خیلی پیش آمده.

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

function get_dir_file_name
hfig = figure('unit','pixel','pos',[200,200,300,100],'menu','none');
uicontrol('pos',[10, 10,70 20],'style','checkbox','string','Directory','tag','dir')
uicontrol('pos',[250, 60,40 30],'string','Go','callback',@getname)
uicontrol('pos',[10, 40,70 20],'style','checkbox','string','Size','tag','size')
uicontrol('pos',[10, 70,70 20],'style','checkbox','string','Dime','tag','date')

function getname(hobject,eventdata,hhandle)
Vdir = get(findall(0,'tag','dir'),'value');
Vsize = get(findall(0,'tag','size'),'value');
Vdate = get(findall(0,'tag','date'),'value');
dir_name = uigetdir;
dir_file = dir(dir_name);
num_file = size(dir_file);
fid=fopen('temp.rtf','w');
controlchar =0;
for i=3:num_file(1)
if ~ dir_file(i).isdir
fprintf(fid,'%s\t',dir_file(i).name);
if Vdate
fprintf(fid,'%s\t',dir_file(i).date);
end
if Vsize
fprintf(fid,'%s\t',num2str(dir_file(i).bytes));
end
fprintf(fid,'%s\n','');
elseif (dir_file(i).isdir & Vdir)
fprintf(fid,'%s\t',dir_file(i).name);
if Vdate
fprintf(fid,'%s\t',dir_file(i).date);
end
fprintf(fid,'%s\t','directory');
fprintf(fid,'%s\n','');
end
end
fclose(fid);
 

 

[ دوشنبه 6 شهریور 1385 - 01:08 ق.ظ ]
[ویرایش شده در : شنبه 21 مهر 1386 - 11:10 ق.ظ]

[ پیام ()|| امین باشی ] [خودم , ] [+]

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

فرض کنید 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;
 

[ یکشنبه 5 شهریور 1385 - 12:08 ب.ظ ]
[ویرایش شده در : یکشنبه 5 شهریور 1385 - 02:08 ق.ظ]

[ پیام ()|| امین باشی ] [حل معادلات دیفرانسیل پاره ای وابسته به زمان در یک بعد , ] [+]

حل عددی معادلات دیفرانسیل با درجات بالاتر از یک

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

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

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

f=dsolve('D2y =(y-6*Dy)/t/4','y(1)=2','y(2)=3');

 ezplot(f,[1 , 10]);

و اما حل عددی :

فرض کنید

 

در نتیجه خواهیم داشت

 

 

و اگر آن را به شکل ماتریس بنویسیم

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

function dy=odefun(t,y)
A = [0 1;1/t/4 -6/t/4];
dy = A*y;

و در خط فرمان مطلب دستور زیر را

[t,y]=ode45('odefun',[1 ,10],[2;3]);

متغییر y دو ستون دارد که ستون اول به y1 و ستون دوم یه y2 اختصاص دارد که y1 جواب معادله و y2 مشتق آن است.

 

 

[ چهارشنبه 1 شهریور 1385 - 11:08 ق.ظ ]
[ویرایش شده در : - - -]

[ پیام ()|| امین باشی ] [حل عددی IVP , ] [+]

+ محاسبه سری فوریه به صورت عددی
+ تعریف هندسه pde
+ معادلات هذلولی ۲
+ معادلات هذلولی
+ رسم pde
+ تبدیل معادلات سمبولیک به عددی
+ حل pde
+ حل تحلیلی معادلات جبری
+ حل معادلات خطی
+ BVP
+ چاپ محتویات یک دایرکتوری
+ pdepe
+ معادلات با درجه بالاتر
+ حل عددی IVP
+ fourier

---

صفحات :

---

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

All right reserved©2005 Amin Bashi