معادله پواسون یک معادله دیفرانسیل جزیی بیضوی (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
% saeed59tb@gmail.com
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)
|
در مورد معادله لاپلاس کافیست قرار دهیم،
حل معادله پواسون دو بعدی به روش ماتریسی
در این مورد با استفاده از روابط (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)
|
|
حل معادله پواسون در دو بعد با استفاده از روش ماتریسی
|
|
|
