Thứ Năm, 20 tháng 2, 2014
Tài liệu CHƯƠNG 9: PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG ppt
407
f=inline(ʹ0ʹ,ʹxʹ,ʹyʹ);
g=inline(ʹ0ʹ,ʹxʹ,ʹyʹ);
x0=0;
xf=4;
Mx=20;
y0=0;
yf=4;
My=20;
bx0=inline(ʹexp(y)‐cos(y)ʹ,ʹyʹ);%(vd.2a)
bxf=inline(ʹexp(y)*cos(4)‐exp(4)*cos(y)ʹ,
ʹyʹ);%(vd.2b)
by0=inline(ʹcos(x)‐exp(x)ʹ,ʹxʹ);%(vd.3a)
byf=inline(ʹexp(4)*cos(x)‐exp(x)*cos(4)ʹ,ʹxʹ);%(vd.3b)
D=[x0xfy0yf];
maxiter=500;
tol=1e‐4;
[U,x,y]=poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,maxiter);
clf
mesh(x,y,U)
axis([0404‐100100])
§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
vànhạnđược:
+
+−
−+ −
=
∆∆
kkkk1k
i1 i i1 i i
2
u2uu u u
A
xt
(2)
408
Côngthứcnàycóthểgóigọnvàothuậttoánsau,gọilàthuậttoánEulẻtiến
tườngminh:
+
+−
∆
=++− =
∆
k1 k k k
ii1i1 i
2
t
ur(uu)(12r)urA
x
(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ó:
ππ
−
π
⎡
⎤
⎞
⎛
λ= + + − = − −
⎜
⎟
⎢
⎥
⎝
⎠
⎣
⎦
j
j
PP
re e (1 2r) 1 2r1 cos
P
(5)
Dotaphảicó|λ|≤1vớibàitoánkhôngcónguồnnênđiềukiệnổnđịnhlà:
∆
=≤
∆
2
t1
rA
x2
(6)
Taxâydựnghàm
fwdeuler()đểthựchiệnthuậttoántrên
function[u,x,t]=fwdeuler(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‐sodoancontheox
%N‐sodiemtheot
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=
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
xt
(7)
−
−+
∆
−++ − = =
∆
kkkk1
i1 i i1 i
2
t
ru (1 2r)u ru u r A
x
(8)
Nếucácgiátrị
k
0
u và
k
M
u ởcảhaiđầuđãchotrướctừđiềukiệnbiênkiểu
Dirichletnênphươngtrình(8)đưatớihệphươngtrình:
−
−
−
−
−−
−
−−
⎡
⎤⎡ ⎤
+
⎡⎤
+−
⎢
⎥⎢ ⎥
⎢⎥
−+−
⎢
⎥⎢ ⎥
⎢⎥
⎢
⎥⎢ ⎥
⎢⎥
−+−
⎢
⎥⎢ ⎥
⎢⎥
=
⎢
⎥⎢ ⎥
⎢⎥
⎢
⎥⎢ ⎥
⎢⎥
+−
⎢
⎥⎢ ⎥
⎢⎥
⎢
⎥⎢ ⎥
⎢⎥
−+
+
⎣⎦
⎣
⎦⎣ ⎦
MMMMMM
MM
L
L
kk1k
110
kk1
22
kk1
33
kk1
M2 M2
kk1k
M1 M1 M
uuru
12rr0000
uu
r12r r 0 0 0
0r12rr00
uu
000 12rr
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)
2x
(10)
vàghépnóvớiphươngtrìnhcóẩn
k
0
u:
−
−
−++ −=
kkkk1
1010
ru (1 2r)u ru u (11)
đểcóđượcphươngtrình:
−
′
+−=− ∆
kkk1
010 0
(1 2r)u 2ru u 2rb (k) x (12)
Kếtquảtacóđượchệphươngtrình:
410
−
−
−
−
−
−
−
⎡
+−
′
⎡
−∆
⎡⎤
⎤
⎢
⎢
⎢⎥
⎥
−+−
⎢
⎢
⎢⎥
⎥
⎢
⎢
⎢⎥
⎥
−+−
⎢
⎢
⎢⎥
⎥
−+
=
⎢
⎢
⎢⎥
⎥
⎢
⎢
⎢⎥
⎥
−
⎢
⎢
⎢⎥
⎥
+−
⎢
⎢
⎢⎥
⎥
⎢
⎢
⎢⎥
⎥
−+
+
⎦
⎣⎦
⎣
⎣
L
L
L
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
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(13
Điểukiệnổnđịnhcủanghiệmlà:
ππ
−
−++−=
λ
jj
PP
1
re (1 2r) re
hay:
λ=
π
⎡⎤
+−
⎢⎥
⎣⎦
1
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
411
A(i,i)=r2;%Pt.(9)
ifi>1
A(i‐1,i)=‐r;
A(i,i‐1)=‐r;end
end
fork=2:N+1
b=[r*u(1,k);zeros(M‐3,1);r*u(M+1,k)]+u(2:M,k‐1);%Pt.(9)
u(2:M,k)=trid(A,b);
end
4.PhươngphápCrank‐Nicholson:Trong(7),xấpxỉđạohàmởvếtráilấyở
thờiđiểmk,trongkhixấpxỉđạohàmởvếphải.Đểcảithiện,talấyđạohàm
ởvếtráilàtrongbìnhcủaxấpxỉđạohàmtại
haiđiểmlàkvàk+1vàcó:
+++ +
+−+−
⎛⎞
−+ −+ −
+=
⎜⎟
∆∆∆
⎝⎠
k1 k1 k1 k k k k1 k
i1 i i1 i1 i i1 i i
22
A u 2u u u 2u u u u
2x x t
(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:
+
+
+
+
−
+
⎡⎤
⎡⎤
+−
⎢⎥
⎢⎥
−+−
⎢⎥
⎢⎥
⎢⎥
⎢⎥
−+−
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
+−
⎢⎥
⎢⎥
⎢⎥
⎢⎥
−+
⎣⎦
⎣⎦
MMMMMM
M
L
L
k1
1
k1
2
k1
3
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
−
⎡⎤
⎡⎤
−
⎢⎥
⎢⎥
−
⎢⎥
⎢⎥
⎢⎥
⎢⎥
−
⎢⎥
⎢⎥
=
⎢⎥
⎢⎥
⎢⎥
⎢⎥
−
⎢⎥
⎢⎥
⎢⎥
⎢⎥
−
⎣⎦
⎣⎦
MMMMMM
M
L
L
k
1
k
2
k
3
k
M1
k
M
u
2(1r)r0000
u
r2(1r)r000
0r2(1r)r00
u
000 2(1r)r
u
000 r2(1r)
u
412
[]
+
⎡⎤
+
⎢⎥
⎢⎥
⎢⎥
⎢⎥
+
⎢⎥
⎢⎥
⎢⎥
′′
⎢⎥
++
⎣⎦
M
k1 k
00
MM
r(u u )
0
0
0
2r b ( k 1) b (k)
(17)
Điềukiệnổnđịnhđượcxácđịnhbằng:
ππ
⎡⎤⎡⎤
⎛⎞ ⎛⎞
λ+ − = − −
⎜⎟ ⎜⎟
⎢⎥⎢⎥
⎝⎠ ⎝⎠
⎣⎦⎣⎦
2 1r1cos 21r1cos
PP
hay:
π
⎡⎤
−−
⎢⎥
⎣⎦
λ= λ≤
π
⎡⎤
+−
⎢⎥
⎣⎦
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);
r2=2*(1+r);
fori=1:M‐1
413
A(i,i)=r2;%Pt.(17)
ifi>1
A(i‐1,i)=‐r;
A(i,i‐1)=‐r;
end
end
fork=2:N+1
b=[r*u(1,k);zeros(M‐3,1);r*u(M+1,k)]
+r*(u(1:M‐1,k‐1)+u(3:M+
1,k‐1))+r1*u(2:M,k‐1);
u(2:M,k)=trid(A,b);%Pt.(17)
end
Đểgiảiphươngtrình:
∂∂
=≤≤≤≤
∂∂
2
2
u(x,t) u(x,t)
0x1,0t0.1
xt
(vd1)
vớiđiềukiệnđầu:
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);
414
figure(3),clf,mesh(t,x,u3)
4.PDEparabolic2chiều:Taxétbàitoánphươngtrìnhviphânđạohàmriêng
parabolichaichiềumôtảsựphânbốnhiệtđộu(x,y,t):
⎡⎤
∂∂ ∂
+=
⎢⎥
∂∂ ∂
⎣⎦
22
22
u(x, y, t) u(x, y, t) u(x, y, t)
A
xy t
(19)
Đểphươngtrìnhcóthểgiảiđượctacầnchođiềukiệnbiên:
=
0
0x
u(x , y,t) b (y,t)
=
f
fx
u(x ,y,t) b (y,t)
=
0
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
xxt
(20)
Taviếtphươngtrìnhtạithờiđiểmtiếptheot
k+1:
+++ ++
+−+−
⎛⎞
−+ −+ −
−=
⎜⎟
⎜⎟
∆∆∆
⎝⎠
k1 k1 k1 k k k k2 k1
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
xxt
(21)
Côngthứcnày,đượcPeacemanvàRachfordđưara,làphươngphápẩnvà
tạonênhệphươngtrình:
(
)
(
)
++ +
−+ −+
− + ++ = − +−
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
∆
=
∆
y
2
t
rA
y
−
∆
f0
x
xx
x=
M
−
∆
f0
y
yy
y=
M
∆
=
T
t
N
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);
ry1=1+2*ry;
ry2=1‐2*ry;
forj=1:Mx‐1%Pt.(22a)
Ay(j,j)=ry1;
ifj>1
Ay(j‐1,j)=‐ry;
Ay(j,j‐1)=‐ry;
end
end
fori=1:My‐1%Pt.(22b)
Ax(i,i)=rx1;
ifi>
1
Ax(i‐1,i)=‐rx;
Ax(i,i‐1)=‐rx;
end
end
416
fork=1:N
u_1=u;
t=k*dt;
fori=1:My+1%Dieukienbien
u(i,1)=feval(bxyt,x(1),y(i),t);
u(i,Mx+1)=feval(bxyt,x(Mx+1),y(i),t);
end
forj=1:Mx+1
u(1,j)=feval(bxyt,x(j),y(1),t);
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:
−
⎡⎤
∂∂ ∂
+=
⎢⎥
∂∂ ∂
⎣⎦
22
4
22
u(x ,y , t) u(x,y , t) u(x ,y , t)
10
xy t
(vd1)
trongmiền:0≤x≤4,0≤y≤4vàtrongkhoảngthơiggian0≤t≤5000
Điềukiệnđầu:
u(x,y,0)=0(vd2a)
vàđiềukiệnbiên:
e
y
cosx‐e
x
cosytạix=0,x=4;y=0,y=4(vd2b)
Đăng ký:
Đăng Nhận xét (Atom)
Không có nhận xét nào:
Đăng nhận xét