Thứ Năm, 20 tháng 2, 2014
Tài liệu CHƯƠNG 7: CÁC PHƯƠNG TRÌNH VI PHÂN THƯỜNG docx
364
h=(xf‐xo)/n;
X=zeros(n+1,1);
M=max(size(yo));%sophuongtrinh(socotcuamatranY)
Y=zeros(n+1,M);
%datdieukiendau
x=xo;
X(1)=x;
y=yo;
Y(1,:)=yʹ;
fori=1:n
ifnargin(fxy)>1
f1
=feval(fxy,x,y);
f2=feval(fxy,x+h,y+h*f1);
else
f1=feval(fxy,x);
f2=feval(fxy,x+h);
end
y=y+h*(f1+f2)/2;
x=x+h;
X(i+1)=x;
Y(i+1,:)=y.ʹ;
end
Đểgiảiphươngtrìnhtadùngchương trình
ctheun.m:
clearall,clc
a=0;
b=1;
y=inline(ʹ2*x+yʹ);
ya=0.5;
n=10;%solantinhchin=10
[x,y]=heun(y,a,b,ya,n)
plot(x,y);
§4.PHƯƠNGPHÁPRUNGE‐KUTTA
365
MặcdùphươngphápHeuntốthơnphươngphápEulernhưngnóvẫn
chưađủđộchínhxácđốivớicácbàitoánthựctế.
XétbàitoánCauchy(1).Giảsửtađãtìmđượcgiátrịgầnđúngy
icủa
y(x
i)vàmuốntínhyi+1củay(xi+1).TrướchếttaviếtcôngthứcTaylor:
)c(y
!m
h
)x(y
!m
h
)x(y
2
h
)x(yh)x(y)x(y
)1m(
1m
i
)m(
m
i
2
ii1i
+
+
+
++⋅⋅⋅+
′′
+
′
+= (11)
vớic∈(x
i,xi+1)và:
[
]
)x(y,x
f
)x(y
iii
=
′
[]
)x(y,xf
dx
d
)x(y
ii
1k
1k
i
)k(
−
−
=
Taviếtlại(11)dướidạng:
)c(y
!m
h
)x(y
!m
h
)x(y
2
h
)x(yhyy
)1m(
1m
i
)m(
m
i
2
ii1i
+
+
+
++⋅⋅⋅+
′′
+
′
=− (12)
TađãkéodàikhaitriểnTaylorđểkếtquảchínhxáchơn.Đểtínhy′
i,y″iv.v.ta
cóthểdùngphươngphápRunge‐Kuttabằngcáchđặt:
)i(
44
)i(
33
)i(
22
)i(
11i1i
krkrkrkryy
+
++=−
+
(13)
trongđó:
⎪
⎪
⎩
⎪
⎪
⎨
⎧
γ+β++=
α++=
=
)kky,bhx(hfk
)ky,ahx(hfk
)y,x(hfk
)i(
2
)i(
1ii
)i(
3
)i(
1ii
)i(
2
ii
)i(
1
(14)
vàtacầnxácđịnhcáchệsốa,b, ;α,β,γ, ;r
1,r2, saochovếphảicủa(13)
khácvớivếphảicủa(12)mộtvôcùngbécấpcaonhấtcóthểcóđốivớih.
KhidùngcôngthứcRunge‐Kuttabậchaitacó:
⎪
⎩
⎪
⎨
⎧
α++=
=
)ky,ahx(hfk
)y,x(hfk
)i(
1ii
)i(
2
ii
)i(
1
(15)
và
)i(
22
)i(
11i1i
krkryy +=−
+
(16)
Tacó:
y′(x)=f[x,y(x)]
[
]
[
]
)x(y,x
f
)x(y,x
f
)x(y
yx
′
+
′
=
′′
Dođóvếphảicủa(12)là:
[]
⋅⋅⋅+
′′
+
′
+ )x(y)y,x(f)y,x(f
2
h
)y,x(hf
iiyiix
2
ii
(17)
Mặtkháctheo(15)vàtheocôngthứcTaylortacó:
366
iii
)i(
1
yh)y,x(h
f
k
′
==
])y,x(
f
k)y,x(
f
ah)y,x(
f
[hk
iiy
)i(
1iixii
)i(
2
⋅
⋅
⋅
+
′
α
+
′
+=
Dođóvếphảicủa(16)là:
⋅
⋅
⋅
+
′
′
α
+
′
++ )]y,x(
f
yr)y,x(
f
ar[h)y,x(
f
)rr(h
iiyi2iix2
2
ii21
(18)
Bâygiờcho(17)và(18)khácnhaumộtvôcùngbécấpO(h
3
)tatìmđượccác
hệsốchưabiếtkhicânbằngcácsốhạngchứahvàchứah
2
:
r
1+r2=1
a.r
1=1/2
α.r
2=1
Nhưvậy: α=a,r
1=(2a‐1)/2a,r2=1/2avớiađượcchọnbấtkì.
Nếua=1/2thìr
1=0vàr2=1.LúcnàytanhậnđượccôngthứcEuler.Nếu
a=1thìr
1=1/2vàr2=1/2.LúcnàytanhậnđượccôngthứcEulercảitiến.
MộtcáchtươngtựchúngtanhậnđượccôngthứcRunge‐Kuttabậc4.
Côngthứcnàyhayđượcdùngtrongtínhtoánthựctế:
k1=h.f(xi,yi)
k
2=h.f(xi+h/2,yi+k1/2)
k
3=h.f(xi+h/2,yi+k2/2)
k
4=h.f(xi+h,yi+k3)
yi
+1=yi+(k1+2k2+2k3+k4)/6
Taxâydựnghàm
rungekutta()đểthựchiệncôngthứcRunge‐Kuttabậc4:
function[x,y]=rungekutta(f,a,b,y0,n)
%Phuong phap Runge‐Kutta de giai phuong trinh yʹ(x) = f(x,y(x)) hay y’ =
%f(x)
ifnargin<4|n<=0
n=100;
end
ifnargin<3
y0=0;
end
y(1,:)=y0(:)ʹ;%
h=(b‐a)/n;
x=a+[0:n]ʹ*h;
ifnargin(f)>1
fork=1:n
367
f1=h*feval(f,x(k),y(k,:));
f1=f1(:)ʹ;
f2=h*feval(f,x(k)+h/2,y(k,:)+f1/2);
f2=f2(:)ʹ;
f3=h*feval(f,x(k)+h/2,y(k,:)+f2/2);
f3=f3(:)ʹ;
f4=h*feval(f,x(k)+h,y(k,:)+
f3);
f4=f4(:)ʹ;
y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6;
end
else
fork=1:n
f1=h*feval(f,x(k));
f1=f1(:)ʹ;
f2=h*feval(f,x(k)+h/2);
f2=f2(:)ʹ;
f3=h*feval(f,x(k)
+h/2);
f3=f3(:)ʹ;
f4=h*feval(f,x(k)+h);
f4=f4(:)ʹ;
y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6;
end
end
Đểgiảiphươngtrìnhtadùngchươngtrình
ctrungekutta.m:
clearall,clc
a=0;
b=1;
y=inline(ʹx+yʹ);
ya=0.5;
n=10;%solantinhchin=10
[x,y]=rungekutta(y,a,b,ya,n)
plot(x,y);
§5.PHƯƠNGPHÁPRUNGE‐KUTTATHÍCHNGHI
368
Vấnđềxácđịnhbướctínhhlà rấtquantrọng.Nếumuốncóđộchính
xáccaothìbướctínhhphảinhỏ.Tuynhiênkhihnhỏ,talạitốnthờigiantính
toán. Hơn nữ
a bước hằng số sẽ không thích hợp trên toàn bộ miền tìm
nghiệm.Vídụnếuđườngcongnghiệmbanđầuthayđổinhanhrồisauđó
gầnnhưkhôngđổithìtaphảidùnghnhỏởđoạnđầu
vàhlớnởđoạnsau.
đâylàchỗmàcácphươngphápthíchnghichiếmưuthế.Chúngđánhgiásai
sốlàmtròntạimỗilầntíchphânvàtựđộnghiệuchỉnhđộlớncủahđểsai
số
nằmtronggiớihạnchophép.
Phương pháp Runge‐Kutta thích nghi còn gọi là phương pháp tích
phânkếthợp.Cáccôngthức nàyđithànhcặp:mộtcôngthứctíchphânbậcm
và
mộtcôngthứctíchphânbậcm+1.Ýtưởnglàdùnghaicôngthứcnàycải
thiệnnghiệmtrongđoạn[x,x+h].Gọikếtquảlày
m(x+h)vàym+1(x+h)tacósai
sốđốivớicôngthứcbậcmlà:
E(h)=y
m+1(x+h)‐ym(x+h)(1)
Chúngtadùngcôngthứckếthợpbậc4và5màđạohàmđượctínhbằng
côngthứcFehlenberg.DovậycôngthứcRunge‐Kuttathíchnghicònđược
gọilàcôngthứcRunge‐Kutta‐Fehlenberg:
=
1
KhF(x,y)
−
=
⎛⎞
=+ +
⎜⎟
⎝⎠
∑
i1
iii,jj
j0
KhFxAh,y BK i=1,2, ,6(2)
=
+= +
∑
6
5ii
i1
y(x h) y(x) CK
(côngthứcbậc5)(3)
=
+= +
∑
6
4ii
i1
y(x h) y(x) DK
(côngthứcbậc4)(4)
Cáchệsốxuấthiệntrongcáccôngthứcnàykhôngduy nhất.Bảngsaucho
cáchệsốtínhtheoCashvàKarp:
i A
i Bi,j Ci Di
1
∗ ∗ ∗ ∗ ∗ ∗
37
378
2825
27648
2
1
5
1
5
∗ ∗ ∗ ∗
0 0
3
3
10
3
40
9
40
∗ ∗ ∗
250
621
18575
48384
4
3
5
3
10
−
9
10
6
5
∗ ∗
125
594
13525
55296
369
5 1
−
11
54
5
2
−
70
27
35
27
∗
0
277
14336
6
7
8
1631
55296
175
512
575
13824
44275
110592
253
4096
512
1771
1
4
Saisốsẽlà:
E(h)=y
5(x+h)‐y4(x+h)=
=
−
∑
6
iii
i1
(C D )K
(5)
ChúýlàE(h)làmộtvectơ,thànhphầnE
i(h)biểudiễnsaisốcủabiếnyi.Sai
sốe(h)tacầnkiểmsoátlà:
=
i
e(h) max E (h) (6)
Tacũngcóthểkiểmsoátsaisốtrungbìnhbìnhphương:
=
=
∑
n
2
i
i1
1
e(h) E ( h)
n
(7)
vớinlàsốphươngtrìnhbậc1.
Việckiểmsoátsaisốđạtđượcbằngcáchthayđổihsaochosai sốtạimỗi
bướctínhậiphỉcỡsaisốmongmuốnε.Saisố
khithựchiênthuậttáonRunge
‐KuttabậcbốnlàO(h
5
)nên:
⎛⎞
≈
⎜⎟
⎝⎠
5
11
22
e(h ) h
e(h ) h
(8)
Giảsửlàtađãtính nghiệmtạibướctínhvớih
1vàcósaisốlàe(h1).Tạibước
tínhvớih
2tamuốncóe(h2)=εthì:
⎡⎤
ε
=
⎢⎥
⎣⎦
1/ 5
21
1
hh
e(h )
(9)
Đểdựphòng,talấy:
⎡⎤
ε
=
⎢⎥
⎣⎦
1/ 5
21
1
h0.9h
e(h )
(10)
Taxâydựnghàm
adaptrk()đểthựchiệnthuậttoánnày:
function[xsol,ysol]=adaptrk(f,xo,x1,y,n)
%TichphanRunge‐Kuttabac5dunggiapphuongtrinhy’=f(x,y)hayy’=
%f(x).
%xo,x1‐doantimnghiem.
%ygiatridau,ndungtimhbandau
370
h=(x1‐xo)/n;
ifsize(y,1)>1;
y=yʹ;%yphailavectohang
end
eTol=1.0e‐9;
n=length(y);
A=[01/53/103/517/8];
B=[00000
1/50000
3/409/40000
3/10‐9/106/500
‐11/545/2‐70/2735/270
1631/55296175/512575/1382444275/110592253/4096];
C=[37/3780250/621125/5940512/1771];
D=[2825/27648018575/4838413525/55296277/143361/4];
%nghiembandau
xsol=zeros(2,1);
ysol=zeros(2,n);
xsol(1)
=xo;
ysol(1,:)=y;
stopper=0;
k=1;
forp=2:5000
%TinhKtu(2)
K=zeros(6,n);
ifnargin(f)>1
K(1,:)=h*feval(f,xo,y);
else
K(1,:)=h*feval(f,xo);
end
fori=2:6
BK=
zeros(1,n);
forj=1:i‐1
BK=BK+B(i,j)*K(j,:);
end
ifnargin(f)>1
371
K(i,:)=h*feval(f,xo+A(i)*h,y+BK);
else
K(i,:)=h*feval(f,xo+A(i)*h);
end
end
%tinhsuthaydoicuaytheo(3)&(4)
dy=zeros(1,n);
E=zeros(1,n);
fori=1:6
dy=dy
+C(i)*K(i,:);
E=E+(C(i)‐D(i))*K(i,:);
end
e=sqrt(sum(E.*E)/n);
%neusaisodatdengiatrichophep,chapnhanketqua
%kiemtrdieukienketthuc
ife<=eTol
y=y+dy;
xo=xo+h;
k=k+1;
xsol(k)=xo;
ysol(k,:)=y;
ifstopper==1;
break
end
end
%tinhlaihtheo(10)
ife~=0;
hnext=0.9*h*(eTol/e)^0.2;
else;
hnext=h;
end
if(h>0)==(xo+hnext>=x1)
hnext=x1‐xo;
stopper=1;
end
h=hnext;
372
end
Đểtìmnghiệmcủaphươngtrìnhviphântadùngchươngtrình
ctadaptrk.m:
clearall,clc
a=0;
b=1;
y=inline(ʹx+yʹ);
ya=0.5;
n=10;%solantinhchin=10
%y=@f4;
[u,v]=adaptrk(y,a,b,ya,n)
plot(u,v)
§6.PHƯƠNGPHÁPBURLIRSCH‐STÖR
1. Phương phápđiểm giữa
: Công thứcđiểm giữa của tích phân số của
′
=yf(x,y)là:
[
]
+= −+y(x h) y(x h) 2hf x,y(x) (1)
Đây là công thức bậc 2, giống như công thức
Euler.Taxemxétphươngphápnàyvìđâylàcơ
sở của phương pháp Burlisch‐Stör dùng tìm
nghiệmcóđộchínhxáccao.Hìnhbênminhhoạ
công thứcđ
iểm giữađối với phương trìnhđơn
dạng
′
=yf(x,y)
.Sựthayđổiytrênhaiphíađược
xácđịnhbằng:
+
−
′
+− −=
∫
xh
xh
y(x h) y(x h) y (x)dx
vàbằngdiệntíchbêndướiđườngcong.Xấpxỉđiểmgiữacủadiện tíchnàylà
diệntíchcủahìnhchữnhậtcógạchchéo.
Bâygiờtaxétưuđiểmcủaphương
phápđiểmgiữakhitìmnghiệmcủa
phươngtrình
′
=yf(x,y)từx=x
0đếnx=x0+Hvớicôngthứcđiểmgiữa.
Tachia đoạntíchphânthànhnđoạnnhỏcóđộdàimỗiđoạnlà
=hH/n
như
hìnhbênvàtính:
=+
10 0
yyhf
=+
20 1
yy2hf
x‐h
h
x+h
x
y’(x)
373
=+
31 2
yy2hf (2)
M
−−
=+
nn2 n1
yy 2hf
Tađãdùngkháiniệmy
i=y(xi)vàfi=f(xi,yi).Ph ương trìnhđầutiêntrong(2)
dùngcôngthứcEulerđểthaychophươngphápđiểmgiữa.Cácphươngtrình
kháclàcáccôngthứcđiểmgiữa.Kếtquảcuốicùnglàtrungbìnhcộngcủay
n
trong(2)vàtacó:
[]
−
+= + +
onn1n
y(x H) 0.5 y (y hf ) (3)
2.NgoạisuyRichardson:Tacóthểthấysaisốtrong(3)là:
=+++L
246
123
Echchch
ĐểgiảmbớtsaisốtadùngphươngphápngoạisuyRichardson.Cụthểtatính
y(x
o+H)vớimộtgiátrịnàođócủahvàrồilặplạiquátrìnhtínhvớih/2.Gọi
kếtquảlàg(h)vàg(h
1)tacóngoạisuyRichardson:
−
+=
1
o
4g(h ) g(h)
y(x H)
3
Taxâydựnghàm
midpoint()đểkếthợpphươngphápđiểmgiữavàphương
phápngoạisuyRichardson.Đầutiênphươngphápđiểmgiữađượcdùngcho
2tíchphân.Sốbướctính đượctănggấpđôitrongcáclầnlặp
sau,mỗilầnlặp
đềudùngngoạisuyRichardson.Chươngtrìnhdừngkhisaisốnhỏh ơnsaisố
chophép.
functiony=midpoint(f,x,x1,y,tol)
%Phuongphapdiemgiuadungchophuongtrinhyʹ=f(x,y)hayy’=f(x).
ifsize(y,1)>1;
y=yʹ;
end%yphailavectohang
ifnargin<5
tol=1.0e‐6;
end
kmax=51;
n=length(y);
r=zeros(kmax,n);
%Batdaubang2buoctichphan
nsteps=2;
r(1,1:n)=mid(f,x,x1,y,nsteps);
rold=r(1,1:n);
xo
x1
x2
x3
xn‐1 xn
H
h
Đăng ký:
Đăng Nhận xét (Atom)
Không có nhận xét nào:
Đăng nhận xét