آموزش مطلب

محاسبه سری فوریه به صورت عددی

امین باشی — 2009-05-28T08:09:24+01:00

قبلا درمورد محاسبه سری فوریه به کمک جعبه ابزار تحلیلی (symbolic) توضیحاتی را داده بودم،امروز نحوه محاسبه سری فوریه به صورت عددی را توضیح می‌دهم.

سری فوریه به شکل‌های گوناگون قابل تعریف است، برای راحتی کار فرم ساده زیر را در بازه‌ی 0< t < 1 را در نظر بگیرید

و ضرایب سری فوریه را می توان به این صورت تعریف کرد

اگر در این قسمت مشکلی دارید به کتابهای ریاضی مهندسی مراجعه کنید.

فرض کنید

برای محاسبه‌ی سری فوریه تنها باید ضرایب آن را حساب کرد، برای محاسبه هرکدام از ضرایب یابدهرکدام از انتگرال‌های فوق را حساب کنیم. برای نمونه من مقدار a_n را حساب می کنم

functiony = anfunc(t)

y = 2*sin(2*pi*t).*sin(2*n*pi*t);

همانطور می دانید هرکدام از این ضرایب یک دنباله از اعدادهستند و با توجه به مقدار n مقدار آنها تغییر می کند. برای محاسبه‌ی تغییرات anfuncبا n 2 راه وجود دارد (مبحث مربوط به ParameterizingFunctions Called by Function Functions را در راهنمای متلبمطالعه کنید). من از تابع تودرتو (nested function) استفاده می‌کنم.

functionmyfunc

forn = 1:10

an(n) = quad(@anfunc,0,1);

end

function y = anfunc(t)

y = 2*sin(2*pi*t).*sin(2*n*pi*t);

end

در ادامهکافی است که ابتدا دیگر ضرایب را حساب کردهو سپس مقدار سری را از جمع کردن آنها محاسبه کنیم.

تعریف هندسه pde

امین باشی — 2007-07-23T00:07:00+01:00

این مقاله حاوی اصطلاحات عجیب و غریب است و خواندن آن برای افراد زیر 7 سال توصیه نمی‌شود

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

http://www.mathworks.com/access/helpdesk/help/toolbox/pde/ug/bqivs1t-1.html

اگر ارل بالا را باز کنید می‌بینید که برای حل pde ابتدا باید ماتریس هندسه مسئله (Geometry Description matrix) را بسازیم. بعد با استفاده از دستور decsg این ماتریس را تجزیه می‌کنیم که به ماتریس حاصل ماتریس هندسه تجزیه شده می‌گویند(Decomposed Geometry matrix). این ماتریس ناحیهای (domain) که معادله بر روی آن حل می شود را مشخص میکند

مثال زیر روش بدست آوردن این دو ماتریس را نشان می‌دهد.

مثال)

s = tf('s');

sys = 1/(s+1)^2;

nyquist(sys)

می بینید که نایکویست این تابع شبیه قلب است؛ و همانطور که می‌دانید نمودار نایکویست در صفحه اعداد مختلط رسم می‌شود. برای حل یک معادله پاره ای را بر روی این شکل باید آن را از فضای اعداد مختلط به فضای اعداد حقیقی تبدیل کنیم(جزییات کار به عهده خودتان).

w=-10:.1:10;

x=(1-w.^2)./(1+w.^2).^2;

y=2*w./(1+w.^2).^2;

plot(x,y)

بد نشد ولی یک مقدار تمیز کاری می‌خواهد.

w = [linspace(-5,-2,7) linspace(-1.9,-.2,20) -.15:.05:.15...

fliplr(-linspace(-1.9,-.2,20)) fliplr(-linspace(-5,-2,7))];

x = (1-w.^2)./(1+w.^2).^2;

y = 2*w./(1+w.^2).^2;

حالا یک مقدار چرخش

xy = [0 x ;0 y ];

XY = [0 1;-1 0]*xy;

اگر به معادلات x,y دقت کنید می بینید هنگامی که w به سمت بینهایت میرود x,y برابر با صفر می شوند ولی در عمل نمی توان w را به سمت بی نهایت میل داد٬ پس صفر را خودمان اضافه کردیم.

هنوز یک مقدر پهن است.

XY = [0.7 0;0 1]*XY;

الان باید ماتریس هندسه را بنویسم

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

در این مثال یک شکل داریم و ماتریس یک ستونی است.

اگر هندسه :

به شکل دایره باشد در سطر اول عدد ۱ ٬ در سطر دوم و سوم به ترتیب مختصات x,y مرکز دایره و در سطر ۴ شعاع دایره را می نویسیم.

به شکل چند ضلعی باشد در سطر اول عدد ۲ ٬ در سطر دوم تعداد اضلاع(n) ٬ در n سطر بعدی xها و در n سطر بعد از آن y ها می نویسیم.

به شکل چهار ضلعی باشد در سطر اول عدد ۳ ٬ را می نویسیم. سطرهای بعدی را مانند حلت چند ظلعی پر می کنیم.

به شکل بیضی باشد در سطر اول عدد ۴ ٬ در سطر دوم و سوم به ترتیب مختصات x,y مرکز بیضی در سطر چهار و پنج اندازه قطر بزرگ و کوچک بیضی و در سطر ششم زاویه چرخش قرار می گیرد.

حالا می‌توانیم ماتریس هندسه بنویسیم.

شکل ما از نوع چند ضلعی است پس در سطر اول عدد 2 را می‌نویسیم.توابع x,y دو تا صفر دارند، یکی در بی‌نهایت و یکی در منفی بی‌نهایت ولی ما فقط یکی را اضافه کردیم چون هر نقطه را فقط 1 بار باید بیاوریم.

x = XY(1,:)';

y = XY(2,:)';

n = length(x);

gd = [2;n;x;y];

dl = decsg(gd);

تا اینجا توانستیم که ماتریس هندسه و ماتریس هندسه تجزیه شده را بدست آوریم. اگر به ارلی در اول مقاله آوردم نگاه کنید می‌بیند که 3 مرحله اول حل الگوریتم معادلات پاره‌ای را انجام داده‌ایم. این 3 مرحله را می‌توان در یک ام-فایل خلاصه کرد. برای نوشتن این ام-فایل از دستور wgeom استفاده می‌کنیم.

fid=wgeom(dl,'cardiodgeom')

این دستور یک ام-فایل با اسم cardiodgeom ایجاد می‌کند. اگر به هر دلیلی این فایل ایجاد نشود مقدار fid برابر با -1 می‌شود.

ادامه دارد...

معادلات هذلولی ۲

امین باشی — 2007-01-28T21:01:00+01:00

مثال مطلب قبلی در نظر بگیرید٬ یک بار دیگر می خواهیم این مثال را حل کنیم.

اولین تغییری که در حل مثال می دهیم٬ نحوه شبکه بندی مسئله است.به نظر می آید که تعداد گره ها برای حل مسئله کافی نباشد با استفاده از دستور refinemesh تعداد گره ها بیشتر می کنیم.با زیاد شدن گره ها تعتداد محاسبات و زمان آن بالا تر می رود ولی ممکن است که نتیجه کار چندان تغییر نکند و یا بهتر باشد که نوع شبکه را عوض کنیم. برای درک بهتر می توانید به کتاب هایی در مورد روش تفاضل محدود-متلب از این روش استفاده می کند- مراجعه کنید.

[p,e,t]=initmesh(g);

[p,e,t]=refinemesh(g,p,e,t);

x=p(1,';

y=p(2,';

u0=atan(cos(pi/2*x));

ut0=3*sin(pi*x).*exp(cos(pi*y));

tlist=linspace(0,5,50);

uu=hyperbolic(u0,ut0,tlist,b,p,e,t,1,0,0,1);

معادلات هذلولی

امین باشی — 2007-01-25T02:01:00+01:00

با استفاده از رابط گرافیکی pdetool حل معادلات پاره ای بسیار ساده است اما بهتر است كه برای درك بهتر، مسئله را به صورت دستی و با نوشتن کد حل کنیم.

فرض کنید بخواهیم معادله موج زیر را حل کنیم :

معادله موج از نوع معادلات هذلولی است و برای حل آن از دستور hyperbolic استفاده می‌كنیم.

u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)

منظور از

u0 و ut0 شرایط اولیه

tlist بازه زمانی حل مسأله

b شرایط مرزی

p,e,t تقسیم بندی مسأله

c,a,f,d ضرایب معادله

می باشد

نوشتن شرایط مرزی را بعدا توضیح می دهم و لی برای این مسئله از pdetool استفاده می كنیم، یك pdetool باز كنید و یك مربع بكشید، شرایط مرزی را مشخص كنید و در متغیرهای b,g ذخیره كیند.

حل نهایی مسئله به این صورت می‌شود.

[p,e,t]=initmesh(g);

x=p(1,:)';

y=p(2,:)';

u0=atan(cos(pi/2*x));

ut0=3*sin(pi*x).*exp(cos(pi*y));

tlist=linspace(0,5,50);

uu=hyperbolic(u0,ut0,tlist,b,p,e,t,1,0,0,1);

for i=1:50

pdesurf(p,t,uu(:,i))

drawnow

end

شاید شما روش من را نپسندید، اینكه یك چیزایی را توضیح نمی دهم و مطالبم گنگ به نظر می‌آید.

به نظر من اگر واقعا به این چیزها احتیاج دارید باید برایش زحمت بكشید.

رسم pde

امین باشی — 2007-01-04T02:01:00+01:00

در مطالب قبلی مقدمات استفاده از pdetoolbox را مرور کردیم و فرض می کنیم که بلد هستید جواب pde را در متغیری مثل u و تقسیم بندی مسئله (mesh) را در متغیر هاب p,e,t ذخیره کنید.

دستور pdemesh برای رسم تقسیم بندی استفاده می شود.

pdemesh(p,e,t)
pdemesh(p,e,t,u)
h=pdemesh(p,e,t)
h=pdemesh(p,e,t,u)

حالات اول دستور٬ تقسیم بندی سیستم را به صورت دو بعدی رسم می کند.

حالت دوم این کار رابه صورت سه بعدی نمایش می دهد٬ ارتفاع هر گره به متغیر u (دما) بستگی دارد.

اما اگر بخواهیم دمای هر گره را بر روی خود آن نشان دهیم به صورت عمل می کنیم:

pdemesh(p,e,t)

hold on

متغیر p دو سطر دارد که سطر اول مولفه x و سطر دوم مولفه y است.تعداد ستون های p به اندازه تعداد گره ها است.

متغیر u یک ستون و به تعداد گره ها سطر دارد. یعنی هر سطر آن متناظر با یک ستون p است.

حالا متغیر u را بتبدیل به رشته گرده و بر روی شکل قرار می دهیم

uu=num2str(u);

nn=length(u);

for i=1:nn

text(p(1,i),p(2,i),uu(i,:))

drawnow

end

دوستی پرسیدن که در هنگام استفاده از plot با پیغام خطای زیر روبرو می شوند:

Attempt to execute SCRIPT plot as a function

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

تبدیل معادلات سمبولیک به عددی

امین باشی — 2006-12-28T23:12:00+01:00

یکی از دوستان سوال کردن چطور می شود معادلات سمبولیک را به عدد تبدیل کرد و گفتند که در نسخه های قبل از متلب ۷ با دستور numeric ٬ این کار قابل انجام بوده است.

من نمی دانم دستور numeric چه کاری انجام می دهد٬ ولی اگر بخواهیم مقدار عددی یک عبارت سبمبولیک را پیدا کنیم از دستور subs استفاده می کنیم

syms x

y = int(x^2)

subs(y,x,3)

قالب کلی اسن دستور به این صورت است

R = subs(S)

R = subs(S,new)

R = subs(S,old,new)

در مثال بالا از حالت سوم دستور استفاده کردیم.

هر عبارت سمبولیک یک متغیر پیش فرض دارد که این متغیر اولین متغیری است که در خروجی دستور findsym دیده می شود.

syms x y z

f = z+y+x

findsym(f)

subs(f,2)

مثال زیر نحوه استفاده از حالات اول دستور را نشان میدهد.

y = dsolve('Dy=y+a')

a = 1000

C1 = 2

subs(y)

حل pde

امین باشی — 2006-12-19T23:12:00+01:00

حل معادلات دیفراسیل جزیی قسمت مهمی از ریاضیات مهندسی است در متلب جعبه ابزار جداگانه ای برای حل اینگونه معادلات قرار داده شده است .مرور کردن این جعبه ابزار به اندازه نوشتن یک کتاب وقت میگیرد.در فایل زیر توضیحات مختصری در مورد نحوه حل این گونه معادلات را نوشته ام.

http://matlabedu.googlepages.com/heat.pdf

حل تحلیلی معادلات جبری

امین باشی — 2006-09-15T21:09:00+01:00

همه ما می توانیم معادله 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

حل معادلات خطی

امین باشی — 2006-09-11T21:09:00+01:00

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

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

BVP

امین باشی — 2006-09-01T22:09:00+01:00

حل معادلات 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 و نحوه نعریف شرایط مرزی داده شود.

چاپ محتویات یک دایرکتوری

امین باشی — 2006-08-27T18:08:00+01:00

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

برنامه زیر این کار را انجام می دهد و فهرست محتویات یم دایرکنپتوری را در یک فایل 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);

pdepe

امین باشی — 2006-08-27T05:08:00+01:00

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

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

معادلات با درجه بالاتر

امین باشی — 2006-08-23T04:08:00+01:00

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

حتما می دانید که هر معادله دیفرانسیل با درجه 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 مشتق آن است.

حل عددی IVP

امین باشی — 2006-08-20T20:08:00+01:00

حل عددی معادلات دیفرانسیل معمولی (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);

fourier

امین باشی — 2006-08-19T22:08:00+01:00

سری فوریه (2)

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

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

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

y=sin(x);

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

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

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

تابع گاما

Y = gamma(A)

تابع خطا

Y = erf(X)