403
CHƯƠNG 9: PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG
§1.KHÁINIỆMCHUNG
Phươngtrìnhviphânđạohàmriêng(PDE)làmộtlớpcácphươngtrình
viphâncósốbiếnđộclậplớnhơn1.Trongchươngnàytasẽkhảosátcác
phươngtrìnhvi
phânđạohàm riêng cấp 2 với hai biếnđộc lập x vày,có
dạngtổngquát:
⎛⎞
∂∂∂ ∂∂
++=
⎜⎟
∂∂∂∂ ∂∂
⎝⎠
222
22
uuu uu
A(x, y) B(x,y) C(x, y) f x, y, u, ,
xxyy xy
(1)
vớix
o≤x≤xf,yo≤y≤yfvàcácđiềukiệnbiên:
=
oyo
u(x,y ) b (x)
=
fyf
∂∂
∇+= + + =
∂∂
22
2
22
u( x ,y) u( x ,y)
u( x ,y) g(x ,y) g( x ,y)u(x ,y) f( x ,y)
xy
(1)
trênmiền
{
}
=≤≤≤≤
ofo f
D(x,y):xxx,yyyvớiđiềukiệnbiêndạng:
==
==
oyo fyf
oxo fxf
u(x, y ) b (x) u(x, y ) b (x)
u( x ,y) b ( y) u( x ,y) b ( y)
(2)
Phươngtrình(1)đượcgọilàphươngtrìnhPoissonnếug(x,y)=0vàgọilà
phươngtrìnhLaplacenếug(x,y)=0vàf(x,y)=0.Đểdùngphươngphápsai
phântachiamiềnthành
Mxđoạn,mỗiđoạndài∆x=(xf‐xo)/Mxdọctheotrục
xvàthànhM
yđoạn,mỗiđoạndài∆y=(yf‐yo)/Mydọctheotrụcyvàthayđạo
hàmbậc2bằngxấpxỉ3điểm:
yx
vớiu
i,j=u(xj,yi) (3.b)
Nhưvậytạimỗiđiểmbêntrong(x
j,xi)với1≤i≤My‐1và1≤j≤Mx‐ậitnhận
đượcphươngtrìnhsaiphân:
+−+−
−+ −+
+==
∆∆
i,j 1 i,j i,j 1 i 1,j i,j i 1,j
i,j i,j i,j
22
u2uu u2uu
gu f
xy
(4)
Trongđó:
u
i,j=u(xj,yi) fi,j=f(xj,yi) gi,j=g(xj,yi)
Cácphươngtrìnhnàysắpxếplạitheocáchnàođóthànhhệphươngtrìnhvới
(M
y‐1)(Mx‐1)biến
{
}
−−−−−
xxyyx
1,1 1,2 1,M 1 2 ,1 2 ,M 1 M 1,2 M 1,M 1
u ,u , ,u ,u , ,u , ,u , ,u .Để
r
2( x y )
∆∆
=
∆+∆
22
xy
22
xy
r
2( x y )
(6)
Bâygiờtakhảosáttiếpcácdạngđiềukiênbiên.CácbàitoánPDEcó2loại
điềukiệnbiên:điềukiênbiênNeumannvàđiềukiênbiênDirichlet.Điềukiện
biênNeumannmôtả
bằng:
=
∂
′
=
∂
o
o
x
xx
u(x,y)
b
(y)
ur(uu)r(u u)r(guf)
+−
′
⎡⎤
=+− ∆+ ++ −
⎣⎦
0
y i,1 i,1 x i x i 1,0 i 1,0 xy i,0 i,0 i,0
ru u 2b(y)x r(u u )r(gu f)
+−
′
⎡
⎤
=+ ++ −− ∆
⎣
⎦
0
y i ,1 x i 1,0 i 1,0 xy i,0 i ,0 i ,0 x i
2r u r ( u u ) r g u f 2 b (y ) x (9)
Nếuđiềukiênbiêntrênbiêndưới(y=y
o)cũnglàkiểuNeumanntasẽviếtcác
phươngtrìnhtươngtựvớij=1,2, ,M
x‐1:
+−
⎡
⎤
′
u2(ruru)rguf2 2
xy
(11)
ĐiềukiệnbiênDirichletchogiátrịhàmtrênbiênnêncóthểthaytrựctiếpvào
phươngtrình.Tacóthểl ấygiátrịtrungbìnhcủacácgiátrịbiênlàmgiátrị
đầucủau
i,j.Taxâydựnghàmpoisson()đểthựchiệnthuậttoánnày:
function[u,x,y]=poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,maxiter)
%giaia(u_xx+u_yy+g(x,y)u=f(x,y)
%trenmienD=[x0,xf,y0,yf]={(x,y)|x0<=x<=xf,y0<=y<=yf}
%voidieukienbien:
%u(x0,y)=bx0(y),u(xf,y)=bxf(y)
%u(x,y0)=by0(x),u(x,yf)=byf(x)
%Mx‐sodoancontrentrucx
%My‐sodoancontrentrucy
%tol:saisochophep
%maxiter:solanlap
x0=D(1);
xf=D(2);
y0=
D(3);
yf=D(4);
dx=(xf‐x0)/Mx;
x=x0+[0:Mx]*dx;
dy=(yf‐y0)/My;
y=y0+[0:My]ʹ*dy;
Mx1=Mx+1;
My1=My+1;
end
end
ifitr>1&max(max(abs(u‐u0)))<tol
break;
end
u0=u;
end
TagiảiphươngtrìnhLaplace:
∂∂
∇= + =
∂∂
22
2
22
u(x, y) u(x, y)
u(x,y) 0
xy
(vd1)
trongmiền0≤x≤4,0≤y≤4vớiđiềukiệnbiên:
=− = −
y
x4
u(0,y) e cosy u(4,y) e cos4 e cosy
(vd2)
==
x4x
§3.PHƯƠNGTRÌNHPARABOLIC
1. Dạng phương trình
: Một phương trình vi phânđạo hàm riêng dạng
paraboliclàphươngtrìnhmôtảsựphânbốnhiệtđộởđiểmxtạithờiđiểmt
củamộtthanh:
∂∂
=
∂∂
2
2
u(x,t) u(x,t)
A
xt
(1)
Đểphươngtrình cóthểgiảiđượctaphảichođiềukiệnbiênu(0, t)=b
0(t),
=
f
fx
u(x , t) b (t)vàđiềukiệnđầuu(x,0)=i0(x)
2.PhươngphápEulertiếntườngminh:Đểápdụngphươngphápsaiphân
hữu hạn, ta chia miên không gian [0, x
f] thành Mđoạn, mỗiđoạn dài
∆=
f
xx/MvàchiathờigianTthànhNphần,mỗiphầnlà∆t=T/N.Sauđóta
thayđạohàmbậc2ởvếtráivàđạohàmbậcởvếph ảicủa(1)bằngcácxấpxỉ
3điểm
(3)
i=1,2, ,M‐1
Đểtìmđiềukiệnổnđịnhcủathuâttoán,tathaynghiệmthử:
π
=λ
j
i
kk
P
i
ue(4)
vớiPlàsốnguyênkháczerovàophươngtrình(3)vàcó:
ππ
−
π
⎡
⎤
⎞
⎛
λ= + + − = − −
⎜
⎟
⎢
⎥
⎝
⎠
⎣
⎦
fori=1:M+1
u(i,1)=it0(x(i));
end
forn=1:N+1
u([1M+1],n)=[bx0(t(n));bxf(t(n))];
end
r=a*dt/dx/dx
r1=
1‐2*r;
fork=1:N
fori=2:M
u(i,k+1)=r*(u(i+1,k)+u(i‐1,k))+r1*u(i,k);%Pt.(3)
409
end
end
3.PhươngphápEulerlùiẩn:Takhảosátmộtthuậttoánkhácgọilàthuật
toánEulerlùi,ẩnsinhradothaythếlùixấpxỉđạohàmđốivớiđạohàmbậc
1trênvếphảicủa(1):
−
+−
−+ −
=
∆∆
kkkkk1
i1 i i1 i i
2
u2uu uu
A
−−
⎡
⎤⎡ ⎤
+
⎡⎤
+−
⎢
⎥⎢ ⎥
⎢⎥
−+−
⎢
⎥⎢ ⎥
⎢⎥
⎢
⎥⎢ ⎥
⎢⎥
−+−
⎢
⎥⎢ ⎥
⎢⎥
=
⎢
⎥⎢ ⎥
⎢⎥
⎢
⎥⎢ ⎥
⎢⎥
+−
⎢
⎥⎢ ⎥
uu
000 r12r
uuru
(9)
ĐiềukiệnbiênNeumann
=
∂
′
=
∂
0
x0
u
b
(t)
x
đượcđưavàophươngtrìnhbằngcách
xấpxỉ:
−
−
′
=
∆
kk
11
0
uu
b
(k)
−
−
⎡
+−
′
⎡
−∆
⎡⎤
⎤
⎢
⎢
⎢⎥
⎥
−+−
⎢
⎢
⎢⎥
⎥
⎢
⎢
⎢⎥
⎥
−+−
⎢
⎢
⎢⎥
⎥
−+
=
⎢
L
MMMMM
MM
L
L
kk
000
kk1
11
kk1
22
kk1
33
kk1
0M2
kk1k
0M1M
12r r 0 0 0 0
uu2rb(k)x
r12r r 0 0 0
uu
0r12rr 00
uu
00 r12r 00
uu
r0
0000 12rr
uu
0000 r12r
uuru
12r1cos
P
λ
≤ 1(14)
Taxâydựnghàm
backeuler()đểthựchiệnthuậttoánnày:
function[u,x,t]=backeuler(a,xf,T,it0,bx0,bxf,M,N)
%Giaiau_xx=u_tvoi0<=x<=xf,0<=t<=T
%Dieukiendau:u(x,0)=it0(x)
%ieukienbien:u(0,t)=bx0(t),u(xf,t)=bxf(t)
%M‐sokhoangcontren
trucx
%N‐sokhoangtheot
dx=xf/M;
x=[0:M]ʹ*dx;
dt=T/N;
t=[0:N]*dt;
fori=1:M+1
u(i,1)=it0(x(i));
end
forn=1:N+1
u([1M+1],n)=[bx0(t(n));bxf(t(n))];
end
r=
a*dt/dx/dx;
r2=1+2*r;
fori=1:M‐1
(15)
vànhậnđượcphươngphápCrank‐Nicholson:
+++
+−+−
∆
−++ − =+− + =
∆
k1 k1 k1 k k k
i1 i i1 i1 i i1
2
t
ru (1 2r)u ru ru (1 2r)u ru r A
x
(16)
VớiđiềukiệnbiênDirichlettạix
0vàđi ềukiệnbiênNeumanntạixMtacóhệ
phươngtrình:
+
+
+
+
−
+
⎡⎤
⎡⎤
+−
⎢⎥
⎢⎥
k1
M1
k1
M
u
2(1r)r0000
u
r2(1r) r 0 0 0
0r2(1r)r00
u
000 2(1r)r
u
000 r2(1r)
u
−
⎡⎤
⎡⎤
−
⎢⎥
⎢⎥
−
⎢⎥
⎢⎥
⎢⎥
⎢⎥
−
⎢⎥
⎢⎥
r2(1r)r000
0r2(1r)r00
u
000 2(1r)r
u
000 r2(1r)
u
412
[]
+
⎡⎤
+
⎢⎥
⎢⎥
⎢⎥
⎢⎥
+
⎢⎥
⎢⎥
⎢⎥
′′
⎢⎥
++
⎣⎦
M
k1 k
00
⎢⎥
⎣⎦
1r1cos
P
1
1r1cos
P
(18)
Taxâydựnghàm
cranknicholson()đểthựchiệnthuậttoántrên:
function[u,x,t]=cranknicholson(a,xf,T,it0,bx0,bxf,M,N)
%Giaiau_xx=u_tvoi0<=x<=xf,0<=t<=T
%Dieukiendau:u(x,0)=it0(x)
%Dieukienbien:u(0,t)=bx0(t),u(xf,t)=bxf(t)
%M‐sokhoang
contrentrucx
%N‐sokhoangtheot
dx=xf/M;
x=[0:M]ʹ*dx;
dt=T/N;
t=[0:N]*dt;
fori=1:M+1
u(i,1)=it0(x(i));
end
forn=1:N+1
u([1M+1],n)=[bx0(t(n));bxf(t(n))];
end
r=a*dt/dx/dx;
r1=2*(1‐r);
u(x,0)=sinπx u(0,t)=0 u(1,t)=0(vd2)
Nhưvậyvới∆x=x
f/M=1/20và∆t=T/N=1/100tacó:
∆
== =
∆
22
t 0.001
rA 1. 0.4
x0.05
(vd3)
Tadùngchươngtrình
ctheat.mđểtìmnghiệmcủa(vd1):
clearall,clc
a=1;%cacthongsocua(vd1)
it0=inline(ʹsin(pi*x)ʹ,ʹxʹ);%dieukiendau
bx0=inline(ʹ0ʹ);
bxf=inline(ʹ0ʹ);%dieukienbien
xf=1;
M=25;
T=0.1;
N=100;
[u1,x,t]=fwdeuler(a,xf,
T,it0,bx0,bxf,M,N);
figure(1),clf,mesh(t,x,u1)
[u2,x,t]=backeuler(a,xf,T,it0,bx0,bxf,M,N);
figure(2),clf,mesh(t,x,u2)
[u3,x,t]=cranknicholson(a,xf,T,it0,bx0,bxf,M,N);
0y
u(x, y , t) b (x, t)
=
f
fy
u(x, y , t) b (x, t)
vàđiềukiệnđầuu(x,y,0)=i
0(x,y)
Tathayđạohàmbậc1theotởvếphảibằngsaiphân3điểmtạiđiểmgiữa
(t
k+1+tk)/2nhưphươngphápCrank‐Nicholson.Tacũngthaythếmộttrong
cácđạohàmbậchaiu
xxvàuyybằngxấpxỉ3điểmtạithờiđiểmtkvàđạohàm
kiatạit
k+1vàcó:
+
+−+−
⎛⎞
−+ −+ −
−=
⎜⎟
⎜⎟
∆∆∆
⎝⎠
kkkkkk k1k
i,j 1 i,j i,j 1 i,j 1 i,j i,j 1 i,j i,j
22
u2uu u2uu uu
A
− + ++ = − +−
k1 k1 k1 k k k
y i 1,j i 1,j y i , j x i , j 1 i ,j 1 x i , j
ru u (1 2r)u ru u (1 2r)u (22a)
với0≤j≤M
x‐1
(
)
(
)
++ + ++ +
−+ −+
− + ++ = − +−
k2 k2 k2 k1 k1 k1
x i ,j 1 i, j 1 x i ,j y i 1,j i 1,j y i ,j
ru u (1 2r)u ru u (1 2r)u (22b)
với0≤i≤M
y‐1
và:
∆
=
∆
x
2
t
rA
x
∆
Taxâydựnghàm
heat2D()đểthựchiệnthuậttoánnày:
function[u,x,y,t]=heat2D(a,D,T,ixy0,bxyt,Mx,My,N)
%Giaiau _t =c(u_xx+u_yy)voiD(1)<=x<=D(2),D(3)<=y<=D(4),0<=t
%<=T
415
%Dieukiendau:u(x,y,0)=ixy0(x,y)
%Dieukienbien:u(x,y,t)=bxyt(x,y,t)voi(x,y)cB
%Mx/My‐cacdoancodoctheotrucx/y
%N‐cackhoangthoigian
dx=(D(2)‐D(1))/Mx;
x=D(1)+[0:Mx]*dx;
dy=(D(4)
‐D(3))/My;
y=D(3)+[0:My]ʹ*dy;
dt=T/N;
t=[0:N]*dt;
%Khoigan
forj=1:Mx+1
fori=1:My+1
u(i,j)=ixy0(x(j),y(i));
end
end
rx=a*dt/(dx*dx);
rx1=1+2*rx;
rx2=1‐2*rx;
ry=a*dt/(dy*dy);
u(My+1,j)=feval(bxyt,x(j),y(My+1),t);
end
ifmod(k,2)==0
fori=2:My
jj=2:Mx;
bx=[ry*u(i,1)zeros(1,Mx‐3)ry*u(i,My+1)]
+rx*(u_1(i‐1,jj)+u_1(i+1,jj))+rx2*u_1(i,jj);
u(i,jj)=trid(Ay,bxʹ
)ʹ;%Pt.(22a)
end
else
forj=2:Mx
ii=2:My;
by=[rx*u(1,j);zeros(My‐3,1);rx*u(Mx+1,j)]
+ry*(u_1(ii,j‐1)+u_1(ii,j+1))+ry2*u_1(ii,j);
u(ii,j)=trid(Ax,by);%Pt.(22b)
end
end
end
Taxétphươngtrình:
−
⎡⎤
∂∂ ∂
+=
⎢⎥
∂∂ ∂
⎣⎦
t]=heat2D(a,D,T,it0,bxyt,Mx,My,N);
mesh(x,y,u)
§4.PHƯƠNGTRÌNHHYPERBOLIC
1. Dạng phương trình
: Phương trình truyền sóng một chiều là PDE dạng
hyperbolic:
∂∂
=
∂∂
22
22
u(x,t) u(x,t)
A
xt
(1)
0≤x≤x
f,0≤t≤T
Điềukiệnbiên:
u(0,t)=b
0(t), =
f
fx
u(x , t) b (t)
vàđiềubiên:
u(x,0)=i
0(x),
=
x
x
M
∆=
T
t
N
vàcóđượcphươngphápsaiphântườngminh:
(
)
+−
+−
=++−−
k1 k k k k1
ii1i1 ii
uruu 2(1r)uu(3)
với:
∆
=
∆
2
2
t
rA
x
uu
i(x)
2t
(5)
vàrútra
−1
i
u đểđưavào(3):
(
)
+−
′
⎡
⎤
=++−−− ∆
⎣
⎦
100 01
ii1i1 ii0i
uru u 2(1r)u u2i(x)t
()
+−
′
=++−+∆
100 0
ii1i1 i0i
1
P
Nhưvậy:
∆
≤=≤
∆
′
π
⎛⎞
−
⎜⎟
⎝⎠
2
2
1t
rrA1
x
1cos
P
Taxâydựnghàm
wave()đểthựchiệnthuậttoántrên:
function[u,x,t]=wave(a,xf,T,it0,i1t0,bx0,bxf,M,N)
%giaiau_xx=u_ttvoi0<=x<=xf,0<=t<=T
%dieukiendau:u(x,0)=it0(x),u_t(x,0)=i1t0(x)
%dieukienbien:u(0,t)=bx0(t),u(xf,t)=bxf(t)
%M‐khoangchiatheox
%
22
22
u(x,t) u(x, t)
xt
(vd1)
0≤x≤2,0≤t≤2
Điềukiệnđầuvàđiềukiệnbiên:
u(x,0)=x(1‐x)
=
∂
=
∂
t0
u
0
t
(vd2a)
=u(0,t) 0 u(1,t)=0(vd2b)
Tadùngchươngtrình
ctwave.mđểgiảiphươngtrìnhnày:
clearall,clc
a=1;
it0=inline(ʹx.*(1‐x)ʹ,ʹxʹ);
i1t0=inline(ʹ0ʹ);%(vd2a)
bx0t=inline(ʹ0ʹ);
bxft=inline(ʹ0ʹ);%(vd2b)
xf=1;
M=20;
T=2;
=
0
x
u(0, y,t) b (y,t)
=
f
fx
u(x , y, t) b (y,t)
=
0
y
u(x,0, t) b (x, t)
=
f
fy
u(x, y , t) b (x, t)
vàđiềubiên:
=
0
u(x, y,0) i (x, y)
=
∂
′
=
∂
0
t0
u(x,y)
i(x,y)
t
y
y
y
M
∆
=
T
t
N
vànhậnđiđếnphươngpháptườngminh:
(
)
(
)
+ −
+− +−
=++−−++−
k1 k k k k k k1
i,j x i,j 1 i,j 1 x y i,j y i 1,j i 1,j i,j
uru u 2(1rr)uru u u (10)
với:
∆
=
∆
2
x
2
+− +−
=++−−++−
100 0001
i,j x i,j 1 i,j 1 x y i,j y i 1,j i 1,j i,j
uru u 2(1rr)uru u u (11)
Nhưvậy,taxấpxỉđiềukiệnđầuvềđạohàmbằngsaiphân:
421
−
−
′
=
∆
11
i,j i,j
0ji
uu
i(x,y)
2t
(12)
vàrútra
−1
i,j
u đểđưavào(11):
()()
+− +−
⎡⎤
=+++
%dieukienbien:u(x,y,t)=
bxyt(x,y,t)voi(x,y)trenbien
%Mx/My‐sokhoangchiatrentrucx/y
%N‐sokhoangchaitheot
dx=(D(2)‐D(1))/Mx;
x=D(1)+[0:Mx]*dx;
dy=(D(4)‐D(3))/My;
y=D(3)+[0:My]ʹ*dy;
dt=T/N;t=[0:N]*dt;
%khoigan
u=zeros(My+1,
Mx+1);
ut=zeros(My+1,Mx+1);
forj=2:Mx
fori=2:My
u(i,j)=it0(x(j),y(i));
ut(i,j)=i1t0(x(j),y(i));
end
end
adt2=a*dt*dt;
rx=adt2/(dx*dx);
ry=adt2/(dy*dy);
422
rxy1=1‐rx‐ry;
rxy2=rxy1*2;
u_1=u;
fork=0:N
t=k*dt;
fori=1:My+1%dieukienbien
Taxétphươngtrình:
⎡⎤
∂∂
∂
+=
⎢⎥
∂∂ ∂
⎣⎦
22
2
222
u(x,y,t) u(x,y,t)
1u(x,t)
4x y t
(vd1)
0≤x≤2,0≤y≤2,0≤t≤2
423
Điềukiệnđầuvàđiềukiệnbiên:
u(0,y,t)=0u(2,y,t)=0 u(x,0,t)=0u(x,2,t)=0 (vd2)
π
π
=
y
x
u(x,y,0) 0.1sin sin
22
∂∂
++ =
∂∂
22
22
u(x,y) u(x,y)
g(x ,y)u(x ,y) f( x ,y)
xy
(1)
trongmiềnDbaobởibiênBvàtrênbiêncócácđiềukiện:
u(x,y)=b(x,y) trênB(2)
CácbướcdùngFEMđểgiảiphươngtrìnhgồm:
)ChiamiềnDthànhNsmiềncon{S1,S2, ,SNs}códạnghìnhtamgiác
)MôtảvịtrícủaNnnútvàđánhsốchúngbắtđầutừcácnúttrênbiên:
n=1,2, ,N
bvàcácnútbêntrong:n=Nb+1,Nb+2, ,Nn
)Xácđịnhcáchàmnộisuy,hìnhdạngvàcơsở:
{
}
φ=φ= ∀∈
nn,s s
(x,y) s 1, ,N (x,y) D (3a)
φ=+ +
n,s n,s n,s n,s
(x,y) p (1) p (2)x p (3)y
chomỗimiềnS
s(3b)
[][]
[
]
[
]
=ϕ+ϕ
TT
11 22
cc(4)
Trongđó:
[]
⎡⎤
ϕ=φ φ φ
⎣⎦
L
b
T
112 N
[]
⎡
⎤
=
⎣
⎦
L
b
T
112 N
ccc c(5a)
∑∑
nn
NN
s n n,s n n,s n,s n,s
i1 i1
(x,y) c (x,y) c p(1)p(x)p(3)y (6)
)Đặtcácgiátrịcủahệsốnútbiêntrong[c1]bằngcácgiátrịbiêntương
ứngvớiđiềukiệnbiên
)Xácđịnhtrịsốcủahệsốnútbêntrongtrong[c2]bằngcáchgiảihệ
phươngtrình:
[A
2][c2]=[d](7)
trongđó:
[]
=
⎧⎫
⎛⎞⎛⎞
∂∂ ∂∂
⎪⎪
⎛⎞⎛⎞
=ϕϕ+ϕϕ∆
⎡⎤ ⎡⎤ ⎡⎤ ⎡⎤
⎨⎬
⎜⎟⎜⎟
⎜⎟⎜⎟
⎣⎦ ⎣⎦ ⎣⎦ ⎣⎦
∂∂ ∂∂
⎝⎠⎝⎠
⎝⎠⎝⎠
⎣⎦
L
b
T
1,s 1,s 2,s N ,s
∂
⎡⎤
ϕ=
⎡⎤
⎣⎦
⎣⎦
∂
L
b
T
1,s 1,s 2,s N ,s
p (2) p (2) p ( 2)
x
∂
⎡⎤
ϕ=
⎡⎤
⎣⎦
⎣⎦
∂
L
22,s2,s2,s2,ss
s1
AS
xx yy
425
{}
=
−ϕϕ∆
⎡⎤⎡⎤
⎣⎦⎣⎦
∑
s
N
T
ss 2,s 2,s s
s1
g(x ,y ) S
(9)
++
⎡⎤
ϕ=φ φ φ
⎡⎤
⎣⎦
⎣⎦
L
bb n
T
2,s N 1,s N 2,s N ,s
p (3) p ( 3) p (3)
y
[]
[][]
=
=− − ϕ ∆
∑
s
N
11 ss2,s
s1
dAc f(x,y)S(10)
(x
s,ys)làtrongtâmcủamiềnconSs
FEMdựatrênnguyêntắclànghiệmcủa(1)cóthểnhậnđượcbằngcách
cựctiểuhoáhàm:
⎧
⎡⎤
∂∂
⎡⎤
=+
⎨
⎢⎥
⎢⎥
∂∂
⎣⎦
⎣⎦
[][][][] [][]
}
−ϕϕ+ ϕ
TT T
g(x,y) c c 2f(x,y) c dxdy (12)
Điềukiệnđểhàmnàycựctiểutheo[c]là:
[]
[]
[][]
{
[][]
∂∂ ∂
=ϕϕ+ϕ
∂∂ ∂
∫∫
TT
2
2
R
d
J
cc
dc x x y
[]
[][]
%N(n,1:2):x&ytoadocuanutthun
426
%S(s,1:3):nutthuscuamiencontamgiacthus
Nn=size(N,1);%tongsonut
Ns=size(S,1);%tongsocacnutcuamiencontamgiac
forn=1:Nn
fors=1:Ns
fori=1:3
A(i,1:3)=
[1N(S(s,i),1:2)];
b(i)=(S(s,i)==n);%hamcosothunbang1chionutthun
end
pnt=A\bʹ;
fori=1:3\
p(n,s,i)=pnt(i);
end
end
end
Đểxácđịnhcácvectơhệsố[c]củanghiệm(4)nhờ(7)vàcácđathức
nghiệm φ
s
(x,y)nhờ(6)đốivớimỗimiềncons=1,2, ,Nstaxâydựnghàm
femcoef():
function[U,c]=femcoef(f,g,p,c,N,S,Ni)
%p(i,s,1:3):cachesocuahamcosophi(i)cuamienconthun
%c=[.11.00.]voicacgiatribienva0voicacnutbentrong
end
d=d‐A12(1:Ni,1:Nb)*c(1:Nb)ʹ;%Pt(10)
c(Nb+1:Nn)=A12(1:Ni,Nb+1:Nn)\d;%Pt.(7)
fors=1:Ns
forj=
1:3
U(s,j)=c*p(:,s,j);%Pt.(6)
end
end
TrướckhidùngFEMđểgiảimộtPDEtaxemthửhàmcơsở(hàmhìnhdạng)
φ
n(x,y)đốivớimỗinútn=1,2, ,Nnđượcđịnhnghĩađốivớitátcảcácmiền
conhìnhtamgiácsaochoφ
nbằng1chỉtạinútnvàbằng0 tạicácnútkhác
đượctạobởihàm
fembasisfth()hoạtđộngnhưthếnào.Tasẽvẽhàmhình
dạngcủamiềnđượcchiathành4tamgiácnhưhìnhsau:
Toạđộcủanút: Sốnútmỗimiềncon
N=[‐11; S=[125;
Copyright: Tài liệu đại học ©