آمار کلی

بازدیدکنندگان : 3009221

Who's Online

ما 67 مهمان آنلاین داریم

آخرين ارسالهاي تالار

پرداخت آنلاین بانک ملت

مبلغ تراکنش (ريال):
حل معادله پواسون -روش تفاضل متناهی مشاهده در قالب PDF چاپ فرستادن به ایمیل
نوشته شده توسط saeed   
سه شنبه, 24 فروردین 1389 ساعت 17:55

حل معادله پواسون -روش تفاضل متناهی ( Finite Differencing Method )

 

معادله پواسون یک معادله دیفرانسیل جزیی بیضوی (Elliptic PDE) است و در دو بعد بشکل زیر نوشته می شود،

(1)

 

با استفاده از تقریب مشتق مرکزی، مشتقات مرتبه دو بصورت زیر در می آیند،

 

(2)

(3)

برای سادگی محاسبات فرض می کنیم h=∆x=Δy در اینصورت با جایگذاری رابطه های (2) و (3) در معادله (1) خواهیم داشت،

(4)

این رابطه در هر نقطه (i,j)Φ واقع شده در مشبندی برقرار است. در حالتیکه g(x,y)=0 (یعنی حالتیکه چشمه بار وجود ندارد) معادله پواسون به معادله لاپلاس تقلیل می یابد در اینصورت رابطه (4) برای معادله لاپلاس بشکل زیر بازنویسی می شود،

 

 

(5)

رابطه (5) نشان می دهد که مقدار Φ برای هر نقطه، میانگین مقادیر Φ ها در چهار نقطه مجاور آن است. سلول محاسباتی پنج نقطه ای توصیف شده در رابطه (5) در شکل زیر نشان داده شده است( اعداد واقع شده در هر سلول مقدار ضریب هر سلول را نشان می دهد)،

اعمال روش تفاضل متناهی (FDM) به مجموعه ای از معادلات جبری می شود. برای یافتن پاسخ مجموعه معادلات دو روش تکرار و ماتریسی مرسومترند.

روش تکرار

روشهای تکرار عموما برای حل سیستم بزرگی از معادلات بکار می روند. شیوه کار در اینگونه موارد آنست که یک تقریب اولیه برای محاسبه تقریب ثانویه استفاده می شود. تقریب ثانویه نیز برای محاسبه سومین تقریب بکار می رود و این روند تا رسیدن به یک همگرایی ادامه می یابد. سه روش تکرار مرسوم ژاکوبی (Jacobi) ، گاوس – سیدل (Gauss - Seidel) و فراواهلش متوالی (successive over-relaxation (SOR)) می باشد.

- روش ژاکوبی

در این روش رابطه (4) را بصورت زیر بازنویسی می شود.

(6)

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

- روش گاوس – سیدل

با دقت کمی متوجه میشویم که در محاسبه (i,j)Φ برای تکرار k+1 اطلاعات مربوط به جملات (Φ(i-1,j و  (Φ(i,j-1 برای تکرار k+1 موجود هستند و از آنجایی این جملات نسبت به جملات مشابه تکرار kام دقیقترند میتوان آنها را جایگزین نماییم. در اینصورت رابطه تکرار گاوس – سیدل بشکل زیر بدست می آید،

(7)

- روش SOR

برای اعمال روش SOR نخست جمله باقیمانده، (R(i,j، در گره iو j بصورت اختلاف بین پاسخ های k و k+1امین تکرار تعریف می کنیم،

 

(8)

توجه داشته باشید k+1(i,j)Φ  از الگوریتم گاوس –سیدل (رابطه 8) جایگزین شده است. مقدار باقیمانده که با R (i,j)نشان داده میشود را می توان بعنوان جمله تصحیح در نظر گرفت که برای نزدیک شدن به پاسخ صحیح باید به (i,j)Φ اضافه می شود. با همگرایی بسمت پاسخ صحیح R(i,j) به صفر میل میکند.. برای تسریع آهنگ همگرایی باقیمانده را در پارامتر ω ضرب می کنیم و آنرا با (i,j)Φ مربوط به kامین تکرار جمع می کنیم تا (i,j)Φ در (k+1)امین تکرار بدست آید. بنابراین رابطه تکرار بصورت زیر در خواهد آمد،

(9)

پارامتر ω فاکتور واهلش نامیده می شود و عددی بین 1 و 2 است. مقدار بهینه ω از طریق آزمون و خطا مشخص می گردد.

اسکریپت زیر معادله پواسون در دو بعد را با روش SOR حل می کند. در اینجا پس از تعداد 200 تکرار برنامه متوقف می شود. الیته می توانید با تعریف یک پارامتر دقت، برنامه را به گونه ای بنویسید که به ازای رسدن به دقت معینی متوقف شود. ناحیه مورد بررسی یک صفحه مستطیلی محصور بین خطوط x=0,x=1,y=0,y=1 می باشد. شرایط اولیه مسئله و همچنین چگالی بار بصورت زیر در نظر گرفته شده است،

 

% finite difference methode of 2-Dimentional poisson equation

% using by SOR methode

clc

clear all

%-----------------------

N=100;   % Number of nodes in each dimention

x=linspace(0,1,N);

y=x;

dx=x(2)-x(1);

dy=dx;

[Y X]=meshgrid(y,x);

V=zeros(N);  % initial guess

% boundry conditions---

V(:,1)=1;      % V(y=0)

V(:,end)=1;    % V(y=1)

V(1,:)=0;      % V(x=0)

V(end,:)=0;    % V(x=1)

%------------------------

rho=X;    % charged density

ep=1;     % epsilon

Vnew=V;

w=1;

tic  % start time

for m=1:200    % iteration loop

for i=2:N-1

for j=2:N-1

Vnew(i,j)=(1-w)*V(i,j)+w/4*(V(i+1,j)+Vnew(i-1,j)...

+V(i,j+1)+Vnew(i,j-1)+1/ep*rho(i,j)*dx^2);

end

 

end

V=Vnew;

end

toc % stop time

mesh(X,Y,V)

xlabel('x')

ylabel('y')

 

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

 

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

روش کار مشابه حالت قبل است , توجه به اینکه مسئله تک بعدیست کارمان ساده تر است. در این مورد رابطه تکرار 7 برای روش فراواهلش SOR بصورت زیر خواهد شد،

(10)

به عنوان مثالی برای این مورد معادله لاپلاس تک بعدی (g(x)=0) را با شرایط اولیه زیر ناحیه x=[0 1] حل می نماییم.

V(x=0)=-1

V(x=1)=2

% finite difference methode of 1-Dimentional Laplas equation

% using by SOR methode

clc

clear all

%-----------------------

N=100;   % Number of nodes in each dimention

x=linspace(0,1,N);

dx=x(2)-x(1);

Vold=zeros(1,N);  % initial guess

% boundry conditions---

Vold(1)=-1;      % V(x=0)

Vold(end)=2;    % V(x=1)

%------------------------

Vnew=Vold;

rho=zeros(1,N);    % charged density

ep=1;           %

w=1.5;

for m=1:1000    % iteration loop

for i=2:N-1

Vnew(i)=(1-w)*Vold(i)+w/2*(Vold(i+1)+Vnew(i-1)...

+1/ep*rho(i)*dx^2);

end

Vold=Vnew;

end

plot(x,Vnew)

 

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

روش ماتریسی

در کوتاهترین جمله متلب یعنی هنربکارگیری ماتریسها. در بسیاری از موارد از نظر ریاضی ((جدا از مسئله سخت افزاری) امکان جایگزینی روشهای ماتریسی  با  حلقه های تکرار وجود دارد و بدون شک متلب بهترین محیط برای این کار است. با بکارگیری قواعد ماتریسی و توابع آماده ماتریسی در متلب امکان نوشتن برنامه هایی عموما کوتاهتر وصد البته بسیار سریعتر (در همگرایی و دقت) را خواهیم داشت. در این قسمت می خواهم به حل ماتریسی معادله پواسون در یک بعد بپردازم. در اینصورت همانگونه که خواهیم دید باید دو ماتریس ایجاد شود. یکی ماتریس ضرایب و دیگری ماتریس معلومها. تشکیل ماتریس ضرایب روی کاغذ و نیز نوشتن برنامه ای برای ایجاد آن در متلب نباید برایتان مشکل باشد و در حقیقت کار ساده ای است. اما هدف من نشان دادن این موضوع است که چگونه با استفاده از توابع و قابلیت های موجود در متلب حداکثر استفاده را برای ایجاد آنچه در ذهن داریم نماییم.

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

 

 

 

 

 

 

 

 

 

 

(11)

همانگونه که پیداست ماتریس ضرایب (ماتریس اول از سمت چپ) یک ماتریس مربعی سه قطری است که درایه های روی قطر اصلی آن -2 و درایه های روی دو قطر دیگر 1 می باشد. چنین ماتریسیهایی را در متلب به اسانی می توان با استفاده از تابع diag خلق نمود. اگر تعداد مشها با درنظرگرفتن نقاط مرزی برابر با N با شد آنگاه دستور زیر این ماتریس     (N-2)*(N-2) را خلق می کند،

A=-2*diag(ones(1,N-2))+diag(ones(1,N-3),1)+diag(ones(1,N-3),-1);

جمله اول عناصر روی قطر اصلی، جمله دوم عناصر بالای قطر اصلی و سومین جمله عناصر پایین قطر اصلی را تولید می کند. با این توضیحات می توانیم اسکریپت مورد نظر را برای مثال قبلی بنویسیم،

clc

clear all

% 1-D poisson equation solution matrix method based finite diffrence

% آدرس ایمیل جهت جلوگیری از رباتهای هرزنامه محافظت شده اند، جهت مشاهده آنها شما نیاز به فعال ساختن جاوا اسکریپت دارید

N=100;

x=linspace(0,2,N);

dx=x(2)-x(1);

phi1=-1;

phiN=1;

A=-2*diag(ones(1,N-2))+diag(ones(1,N-3),1)+diag(ones(1,N-3),-1);

gx=x;

B=dx^2*gx(2:end-1)';

B(1)=B(1)-phi1;

B(end)=B(end)-phiN;

phi=A\B;

phi=[phi1 phi' phiN];

plot(x,phi)

در مورد معادله لاپلاس کافیست قرار دهیم،

gx=zeros(1,N);

 

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

در این مورد با استفاده از روابط (2)  و (3) و جایگذاری آنها در معادله پواسون دوبعدی به معادله ماتریسی مشابه رابطه (9) میرسیم. مهمترین مشکل تشکیل ماتریس ضرایب در حالت دو بعدیست که به نسبت حالت یک بعدی کار مشکلییست. اسکریپت زیر معادله پواسون را به روش ماتریسی برای مسئله ای با همان شرایط اولیه که در بالا انرا به روش SOR حل کردیم، حل می نماید.

ایراد روش ماتریس بخصوص در مواردی که با ابعاد بیش از یک بعد سرو کار داریم بحث حافظه مورد نیاز است. در چنین مواردی بهتر است از روشهایی با دقت بالاتر استفاده نماییم تا مجبور نیاشیم ابعاد ماتریسها را خیلی بزرگ نماییم. در مورد معادلاتی نظیر معادله پواسون در ابعاد بیش از یک بعد روش المان محدود جایگزین مناسبی برای روش تفاضل محدود است چرا که  از دقت بالاتری برخوردار است. ایراد دیگر روش تفاضل متناهی در مواجه با مواردی است که با نواحی مرزی با شکلهای پیچیده تر روبرو هستیم که در این موارد این روش عملا کارایی ندارد. در چنین مواردی روش المان را باید بکار برد. در بخشهای بعدی روش استفاده از جعبه ابزار (toolbox) المان محدود را برای حل چنین مسائلی شرح خواهم داد.

clc

clear; % clear workspace to start new

close; % close previous figures to start new

%------------------

Nx=50;       % number of mesh- x dirextion

Ny=50;       % number of mesh- y dirextion

nx=Nx-2;     % number of mesh- x dirextion without boundry point clc

clear; % clear workspace to start new

close; % close previous figures to start new

%------------------

Nx=50;       % number of mesh- x dirextion

Ny=50;       % number of mesh- y dirextion

nx=Nx-2;     % number of mesh- x dirextion without boundry point x=0 & x=Lx

ny=Ny-2;     % number of mesh- y dirextion without boundry point y=0 & x=Ly

x=linspace(0,1,Nx);

y=linspace(0,1,Ny);

dx=x(2)-x(1);

dy=y(2)-y(1);

X=repmat(x(2:end-1),1,ny)';

Y=repmat(y(2:end-1),nx,1);

Y=Y(:);

rho=X;

ep=1;

a=diag(ones(1,nx));

b=diag(ones(1,nx-1),1);

c=diag(ones(1,nx-1),-1);

A1=(-2*a+b+c)/dx^2-2*a/dy^2;

A=A1;

for i=1:ny-1

A=blkdiag(A,A1);

end

A=A+diag(ones(1,nx*(ny-1)),nx)/dy^2+diag(ones(1,nx*(ny-1)),-nx)/dy^2;

phi0=zeros(1,nx*ny)';

phi0(1:nx:end)=1/dx^2;

phi0(nx:nx:end)=1/dx^2;

phi0(1:nx)=0;

phi0(end-nx+1:end)=0;

B=-rho/ep-phi0;

phi=A\B;

phi=reshape(phi,nx,ny);

surf(phi)

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

 

 

نظر ها (33)
  • ناشناس
    دستتون درد نکنه :wink:
  • ناشناس
    kheili mersi
  • gfhgjj
    thanxxxxxxxxxxxxx
  • پدرام  - تقدیر و تشکر
    با تشکر !!!!!!
  • رقیه  - تشکر
    ممنون
  • ناشناس
    خیلی باحالی ممنون دوستت دارم :wink:
  • ناشناس
    اي ول واقعا كارتون درسته
  • علی نیک فال
    سلام
    وب پرمحتوا و پرباری دارید و واقعا خدا قوت
    قابل توجه دانش پژوهان عرصه مخابرات،رایانه و دانشجویان عزیز و شرکت های فعال در زمینه های مختلف رشته مهندسی مخابرات:

    طراحی ، پیاده سازی و شبیه سازی پروژه های درسی، تحقیقاتی و صنعتی با

    نرم افزار OPNET Modeler

    به همرا دوره آموزش کوتاه مدت

    طرح ها و پروژه های انجام شده:

    1-طراحی و پیاده سازی انواع شبکه های LAN، MAN و WAN

    2- شبکه های PSTN

    3-شبکه UMTS

    4-شبکه WiMAX

    5-MPLSو ATM

    6- شبکه NGN

    7-امنیت شبکه های نسل جدید

    و...

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

    ۱-دوره آموزش رایگان به ازای سفارش کار

    ۲-ارایه نرم افزار برای سفارش دهندگان پروژه به صورت رایگان

    ۳- انجام پروژه در کمترین زمان ممکن با توجه به نیاز مشتری به همراه توضیحات کار های انجام یافته

    ۴- ارایه تحلیل های آماری حاصل از شبیه سازی انجام یافته

    ۵- ارایه مشاوره رایگان

    ۶- ارایه کد های منبع برنامه

    زمینه طراحی شبکه شامل موارد زیر می ب...
  • فاطمه
    عالی بود ممنون
  • ناشناس
    برای حل معادلات وبه دست اوردن ریشه به روش نصف کردن وروش ss در متلب دستور اماده وجود دارد یا اینکه باید خودمان دستورش را بنویسیم :?: :?: مثلا y=(x/2)-sinxj,تو بازه(0,2) باشرط E=0.01
    :?:
  • saeed
    از fsolve استفاده کنید. این دستور نزدیکترین ریشه را به مقدار پیش فرض در بازه مورد نظر جستجو میکند.
  • سعید
    مهندس واقعا عالی بود ومن به شخصه برای شما آرزوی پیروزی می کنم :)
تنها کاربران عضو شده می توانند نظر ارسال کنند!
آخرین بروز رسانی در چهارشنبه, 26 آبان 1389 ساعت 21:42