360
CHƯƠNG 7: CÁC PHƯƠNG TRÌNH VI PHÂN THƯỜNG

§1.BÀITOÁNCAUCHY
 Một phương trình vi phân cấp 1 có thể viết dưới dạng giảiđược
′
=yf(x,y)
màtacóthểtìmđượchàmytừđạohàmcủanó.Tồntạivôsố
nghiệmthoảmãnphươngtrìnhtrên.Mỗinghiệm phụthuộcvàomộthằngsố
tuỳý.Khichotrướcgiátrịbanđầ
ucủaylàyotạigiátrịđầuxotanhậnđược
mộtnghiệmriêngcủaphươngtrình.BàitoánCauchy(haybàitoáncóđiều
kiệnđầu)tóm lạinhưsau:choxsaochob
≥x≥a,tìmy(x)thoảmãnđiều
kiện:
⎩
⎨
⎧
α=
=
′
)a(y
)y,x(
f
)x(y
(1)
 Ngườitachứngminhrằngbàitoánnàycómộtnghiệmduynhấtnếuf
thoảmãnđiềukiệnLipschitz:

2121

=
′
)a(Y
)X,x(
f
)x(Y

với:
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎛
′
′
′
=
′
n
f
f
F

⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎛
=
n
2
1
yy

⎨
⎧
=
′
=
′
)v,u,x(gv
vu

vớiđiềukiệnđầu:u(a)=αvàv(a)=β
 Các phương
 pháp giải phương trình vi phânđược trình bày trong
chươngnàylàcácphươngpháprờirạc:đoạn[a,b]đượcchiathànhnđo ạn
nhỏbằngnhauđượcgọilàcácʺbướcʺtíchphânh=(
b‐a)/n.

§2.PHƯƠNGPHÁPEULER
Giảsửtacóphươngtrìnhviphân:

⎩
⎨
⎧
α=
=
′
)a(y
)y,x(
f
)x(y
(1)

2
i1i
ii1ii1i

Nếu(x
i+1‐xi)khábéthìtacóthểbỏquacács ốhạng(xi+1‐xi)
2
vàcácsố
hạngbậccao
y(x
i+1)=y(xi)+(xi+1‐xi)y′(xi)
Trườnghợpcácmốccáchđều:
(x
i‐1‐xi)=h=(x‐xo)/n
thìtanhậnđượccôngthứcEulerđơngiản:
 y
i+1=yi+hf(xi,yi)   (2)
Vềmặthìnhhọctathấy(1)chokếtquảcàng
chínhxácnếubướchcàngnhỏ.
Taxâydựnghàm
euler()đểthựchiệnthuậttoántrên:
y
x
x
i
xi+1
yi
yi+1

362

dy(1)=y(2)*y(3);
dy(2)=‐y(1)*y(3);
dy(3)=‐0.51*y(1)*y(2);


Đểgiảiphươngtrìnhchobởihàmf1(x,y)tadùngchươngtrình
cteuler.m:

clearall,clc
a=0;

363
b=1;
y=@f1;
ya=[011]ʹ;
m=200;
[x,y]=euler(y,a,b,ya,m)
plot(x,y);

§3.PHƯƠNGPHÁPHEUN
 Phương pháp Heun cònđược gọi là phương pháp hình thang hay
phươngpháp.Chophươngtrình:
 y’=f(t,y)
Tacó:
+
+
+
=−=
∫
k1

2

Vếphải(RHS)củaphươngtrìnhnàycóy
k+1làgáitrịchưabiếttạithờiđiểmtk.
Đểgiảiquyếtvấnđềnàytathayy
k+1ởRHSbằngcôngthứcxấpxỉ:

+
≅+
k1 k k k
yyhf(t,y)
Nhưvậy:

[]
{
}
++
=+ + +
k1 k k k k1 k k k
h
y y f(t ,y ) f (t ,y hf(t ,y )
2

ĐâychínhlàcôngthứcHeun.
Taxâydựnghàm
heun()đểthựchiệnthuậttoántrên:

function[X,Y]=heun(fxy,xo,xf,yo,n)
%Giaiphuongtrinhyʹ(x)=f(x,y(x))hayy’=f(x)
%dungthuattoanHeunvoinbuoctinh

Đểgiảiphươngtrìnhtadùngchương trình
ctheun.m:

clearall,clc
a=0;
b=1;
y=inline(ʹ2*x+yʹ);
ya=0.5;
n=10;%solantinhchin=10
[x,y]=heun(y,a,b,ya,n)
plot(x,y);

§4.PHƯƠNGPHÁPRUNGE‐KUTTA

365
MặcdùphươngphápHeuntốthơnphươngphápEulernhưngnóvẫn
chưađủđộchínhxácđốivớicácbàitoánthựctế.
XétbàitoánCauchy(1).Giảsửtađãtìmđượcgiátrịgầnđúngy
icủa
y(x
i)vàmuốntínhyi+1củay(xi+1).TrướchếttaviếtcôngthứcTaylor:
)c(y
!m
h
)x(y
!m
h
)x(y
2
h

)x(y,xf
dx
d
)x(y
ii
1k
1k
i
)k(
−
−
= 
Taviếtlại(11)dướidạng:
)c(y
!m
h
)x(y
!m
h
)x(y
2
h
)x(yhyy
)1m(
1m
i
)m(
m
i
2

⎪
⎪
⎨
⎧
γ+β++=
α++=
=

)kky,bhx(hfk
)ky,ahx(hfk
)y,x(hfk
)i(
2
)i(
1ii
)i(
3
)i(
1ii
)i(
2
ii
)i(
1
(14)
vàtacầnxácđịnhcáchệsốa,b, ;α,β,γ, ;r
1,r2, saochovếphảicủa(13)
khácvớivếphảicủa(12)mộtvôcùngbécấpcaonhấtcóthểcóđốivớih.
KhidùngcôngthứcRunge‐Kuttabậchaitacó:


]
)x(y,x
f
)x(y,x
f
)x(y
yx
′
+
′
=
′′


Dođóvếphảicủa(12)là:
[]
⋅⋅⋅+
′′
+
′
+ )x(y)y,x(f)y,x(f
2
h
)y,x(hf
iiyiix
2
ii
(17)
Mặtkháctheo(15)vàtheocôngthứcTaylortacó:


′
+= 
Dođóvếphảicủa(16)là:

⋅
⋅
⋅
+
′
′
α
+
′
++ )]y,x(
f
yr)y,x(
f
ar[h)y,x(
f
)rr(h
iiyi2iix2
2
ii21
(18)
Bâygiờcho(17)và(18)khácnhaumộtvôcùngbécấpO(h
3
)tatìmđượccác
hệsốchưabiếtkhicânbằngcácsốhạngchứahvàchứah
2
:

ifnargin<4|n<=0
n=100;
end
ifnargin<3
y0=0;
end
y(1,:)=y0(:)ʹ;%

h=(b‐a)/n;
x=a+[0:n]ʹ*h;
ifnargin(f)>1
fork=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
fork=1:n
f1=h*feval(f,x(k));
f1=f1(:)ʹ;
f2=h*feval(f,x(k)+h/2);

nghiệm.Vídụnếuđườngcongnghiệmbanđầuthayđổinhanhrồisauđó
gầnnhưkhôngđổithìtaphảidùnghnhỏởđoạnđầu
vàhlớnởđoạnsau.
đâylàchỗmàcácphươngphápthíchnghichiếmưuthế.Chúngđánhgiásai
sốlàmtròntạimỗilầntíchphânvàtựđộnghiệuchỉnhđộlớncủahđểsai
số
nằmtronggiớihạnchophép.
 Phương pháp Runge‐Kutta thích nghi còn gọi là phương pháp tích
phânkếthợp.Cáccôngthức nàyđithànhcặp:mộtcôngthứctíchphânbậcm
và
mộtcôngthứctíchphânbậcm+1.Ýtưởnglàdùnghaicôngthứcnàycải
thiệnnghiệmtrongđoạn[x,x+h].Gọikếtquảlày
m(x+h)vàym+1(x+h)tacósai
sốđốivớicôngthứcbậcmlà:
 E(h)=y
m+1(x+h)‐ym(x+h)(1)
Chúngtadùngcôngthứckếthợpbậc4và5màđạohàmđượctínhbằng
côngthứcFehlenberg.DovậycôngthứcRunge‐Kuttathíchnghicònđược
gọilàcôngthứcRunge‐Kutta‐Fehlenberg:

=
1
KhF(x,y)


−
=
⎛⎞
=+ +
⎜⎟

∗ ∗ ∗ ∗ ∗ ∗
37
378

2825
27648
2
1
5

1
5

∗ ∗ ∗ ∗
0 0
3
3
10

3
40

9
40

∗ ∗ ∗
250
621

18575


−
70
27

35
27

∗
0
277
14336
6
7
8

1631
55296

175
512

575
13824

44275
110592
253
4096


n
2
i
i1
1
e(h) E ( h)
n
(7)
vớinlàsốphươngtrìnhbậc1.
Việckiểmsoátsaisốđạtđượcbằngcáchthayđổihsaochosai sốtạimỗi
bướctínhậiphỉcỡsaisốmongmuốnε.Saisố
khithựchiênthuậttáonRunge
‐KuttabậcbốnlàO(h
5
)nên:

⎛⎞
≈
⎜⎟
⎝⎠
5
11
22
e(h ) h
e(h ) h
(8)
Giảsửlàtađãtính nghiệmtạibướctínhvớih
1vàcósaisốlàe(h1).Tạibước
tínhvớih
2tamuốncóe(h2)=εthì:

%f(x).
%xo,x1‐doantimnghiem.
%ygiatridau,ndungtimhbandau

370
h=(x1‐xo)/n;
ifsize(y,1)>1;
y=yʹ;%yphailavectohang
end
eTol=1.0e‐9;
n=length(y);
A=[01/53/103/517/8];
B=[00000
1/50000

3/409/40000
3/10‐9/106/500
‐11/545/2‐70/2735/270
1631/55296175/512575/1382444275/110592253/4096];
C=[37/3780250/621125/5940512/1771];
D=[2825/27648018575/4838413525/55296277/143361/4];
%nghiembandau
xsol=zeros(2,1);
ysol=zeros(2,n);
xsol(1)
=xo;
ysol(1,:)=y;
stopper=0;
k=1;
forp=2:5000

%neusaisodatdengiatrichophep,chapnhanketqua
%kiemtrdieukienketthuc
ife<=eTol
y=y+dy;
xo=xo+h;

k=k+1;
xsol(k)=xo;
ysol(k,:)=y;
ifstopper==1;
break
end
end
%tinhlaihtheo(10)
ife~=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ìmnghiệmcủaphươngtrìnhviphântadùngchươngtrình

=yf(x,y)
.Sựthayđổiytrênhaiphíađược
xácđịnhbằng:

+
−
′
+− −=
∫
xh
xh
y(x h) y(x h) y (x)dx
vàbằngdiệntíchbêndướiđườngcong.Xấpxỉđiểmgiữacủadiện tíchnàylà
diệntíchcủahìnhchữnhậtcógạchchéo.
 Bâygiờtaxétưuđiểmcủaphương
phápđiểmgiữakhitìmnghiệmcủa
phươngtrình
′
=yf(x,y)từx=x
0đếnx=x0+Hvớicôngthứcđiểmgiữa.
Tachia đoạntíchphânthànhnđoạnnhỏcóđộdàimỗiđoạnlà
=hH/n
như
hìnhbênvàtính:

=+
10 0
yyhf
 =+
20 1

2.NgoạisuyRichardson:Tacóthểthấysaisốtrong(3)là:

=+++L
246
123
Echchch

ĐểgiảmbớtsaisốtadùngphươngphápngoạisuyRichardson.Cụthểtatính
y(x
o+H)vớimộtgiátrịnàođócủahvàrồilặplạiquátrìnhtínhvớih/2.Gọi
kếtquảlàg(h)vàg(h
1)tacóngoạisuyRichardson:

−
+=
1
o
4g(h ) g(h)
y(x H)
3

Taxâydựnghàm
midpoint()đểkếthợpphươngphápđiểmgiữavàphương
phápngoạisuyRichardson.Đầutiênphươngphápđiểmgiữađượcdùngcho
2tíchphân.Sốbướctính đượctănggấpđôitrongcáclầnlặp
sau,mỗilầnlặp
đềudùngngoạisuyRichardson.Chươngtrìnhdừngkhisaisốnhỏh ơnsaisố
chophép.

functiony=midpoint(f,x,x1,y,tol)

%kiemtrahoitu.
dr=r(1,1:n)‐rold;
e=
sqrt(dot(dr,dr)/n);
ife<tol;y=r(1,1:n);
return;
end
rold=r(1,1:n);
end
error(ʹPhuongphapdiemgiuakhonghoituʹ)

functiony=mid(f,x,xf,y,nsteps)
%Congthucdiemgiua
h=(xf‐x)/nsteps;
y0=y;
ifnargin(f)
>1
y1=y0+h*feval(f,x,y0);
else
y1=y0+h*feval(f,x);
end
fori=1:nsteps‐1
x=x+h;
ifnargin(f)>1
y2=y0+2.0*h*feval(f,x,y1);
else
y2=y0+2.0*h*feval(f,x);
end

y0=y1;

%H=dotangsaumoilantinh
ifsize(y,1)>1
y=yʹ;
end%yphailavectohang
ifnargin<6
tol=1.0e‐8;
end
n=length(y);
xout=zeros(2,1);
yout=zeros(2,n);
xout(1)=x;
yout(1,:)=y;
k=1;
whilex<x1

376
k=k+1;
H=min(H,x1‐x);
y=midpoint(f,x,x+H,y,tol);
x=x+H;
xout(k)=x;
yout(k,:)=y;
end


Đểgiảiphươngtrìnhtadùngchương trình
ctburlischstoer.m:

clearall,clc
a=0;

Ey()hxxh
(m 1)!

Dùngxấpxỉđạohàm:

+
+−
ξ≈
(m) (m)
(m 1)
y(xh)y(x)
y()
h

tacó:

⎡⎤
=+−
⎢⎥
+
⎣⎦
m
(m) (m)
h
Ey(xh)y(x)
(m 1)!
(2)
Taxâydựnghàm
taylor()đểgiảibàitoántrên:


Tadùngchươngtrình
cttaylor.mđểgiảiphươngtrình:

clearall,clc
y=@f5;
a=0;
b=2;
ya=[01];
h=0.2;
[x,y]=taylor(y,a,ya,b,h)

378
plot(x,y);

§8.PHƯƠNGPHÁPDỰĐOÁN‐HIỆUCHỈNH
1.PhươngphápAdam‐Bashfort‐Moulton
:Năm1855,nhàtoánhọcngười
AnhAdamsđềxuấtmộtphươngphápđabướcgiảibàitoánCauchytheoyêu
cầucủaôngBashforth,mộtchuyêngiakỹthuậtpháobinhAnh.Kếtquảcủa
Adams
saunàybịquênlãng.Mãiđếnđầuthếkỷ20,nhàtoánhọcNauykhi
tínhquỹđạocủahạtđiệntíchrờixamặttrờivớivậntốclớnđãphátminhlại
công thức
Adams. Sau này viện sỹ Krylovđã hoàn thi ện phương pháp
Adams. Phương pháp Adams‐Bashfort‐Moulton (ABM) gồm hai bước.
Bướcdầutiênlàxấpxỉf(x,y)bằngmộtđathức(vídụđathứcLagrange)b
ậc4
qua4điểm:

(

)
(
)
(
)
(
)
{
}
−− −− ++k2 k2 k1 k1 k k k1 k1
t,f ,t,f ,t,f,t,f 
vànhậngiátrịhiệuchỉnh:
+++
=
k1 k1 k1
ff(t,p)

(
)
−− +
′
=+ =+ − + +
∫
h
kk 3 k k2 k1 k k1
0
h
c y l (t)dt y f 5f 19f 9f
24
  (1b)

23!
(2b)
vàthaythếcácđạohàmbậc1,2,3bằngcácxấpxỉ

−−−
+
⎛⎞
−+−+
⎜⎟
=+ + + +
⎜⎟
⎜⎟
⎝⎠
L
2
k3 k2 k1 k
3(4)
k1 k k k
13 11
ff3f f
h1
32 6
yyhf hf
2h4


379

−−−
−+ − +

4! h 2 120


(
)
−−−
=+ − + − + + +L
5(4)
kk3k2k1k k
h251
y9f37f59f55fhf
24 720


+
≈+
5(4)
k1 k
251
phf
720
(3a)

−− +
++ +
⎛⎞
−+−+
⎜⎟
=+ − + +
⎜⎟

hf
3! h 12



−− +
++
−+ −+
⎛⎞
−++++
⎜⎟
⎝⎠
LL
45
(4) (4)
k2 k1 k k1
k1 k1
3
h f 3f 3f f 3 h
hf f
4! h 2 120


(
)
−− + +
=+ − + + − +L
5(4)
k k2 k1 k k1 k1
h19

+++ +
=−≈− ≅−
5(4) 5
C,k1 k1 k1 k1
19 19
Eyc hf Kh
720 720
(4b)
DoKchưabiếtnêntaphảitìmnó.Tacó;

++++ + +
−=−≅ ≡ ≡−
5
P,k1 C,k1 k1 k1 P,k1 C,k1
270 270 270
EEcp Kh E E
720 251 19
 (5)
Dovậytacócáccôngthứcdùngđểđánhgiásaisố:

()
+++ ++
=−≅ −
P,k1 k1 k1 k1 k1
251
Eyp cp
720
(6a)

()

()
−− ++
⎡
⎤
=+ − + +
⎢
⎥
⎣
⎦
kk k2 k1 k k1k1
h
c y f 5f 19f 9f t ,m
24
 (7c)

()
++ ++
=− −
k1 k1 k1 k1
19
yc cp
270
(7d)
Taxâydựnghàm
odeabm()đểthựchiệnthuậttoánnày:

function[t,y]=odeabm(f,to,tf,y0,n)
%PhuongphapAdams‐Bashforth‐Moulton
%degiaiptyʹ(t)=f(t,y(t))hayyʹ(t)=f(t)
if(nargin<5)|(n<0)

c1=y(k,:)+h242*[F(2:4,:);feval(f,t(k+1))ʹ];%Pt.(7c)
end
y(k+1,:)=c1‐KC12*(c1‐p1);%Pt.(7d)
p=p1;
c=c1;%capnhatcacgiatridudoan/hieu
chinh
ifnargin(f)>1
F=[F(2:4,:);feval(f,t(k+1),y(k+1,:))ʹ];
else
F=[F(2:4,:);feval(f,t(k+1))ʹ];
end
end


Đểgiảiphươngtrìnhtadùngchươngtrình
ctodeabm.m:

clearall,clc
a=0;
b=1;
y=@f1;
ya=[011]ʹ;
n=10;
[t,y]=odeabm(y,a,b,ya,n)
plot(t,y)


2. Phương pháp Hamming: Thuật toán Hamming cũng như thuật toán
Adams‐Bashforth‐Moultonnhưngcáccôngthứcdựbáo/hiệuchỉnhlà:
Dựđoán:

c0.1259yy 3hf 2fft,m
 (8c)

()
++ + +
=− −
k1 k1 k1 k1
9
yc cp
121
(dd)
Taxâydựnghàm
hamming()đểthựchiệnthuậttoánnày:
function[t,y]=hamming(f,to,tf,y0,n)

382
%PhuongphapHammingdegiaiphuongtrinhyʹ(t)=f(t,y(t))hayyʹ=f(t)
if(nargin<5)|(n<=0)
n=100;
end
h=(tf‐to)/n;
ts=to+3*h;
[t,y]=rungekutta(f,to,ts,y0,3);%KhoiganbangRunge‐Kutta
t=
[t(1:3)ʹt(4):h:tf];
fork=2:4
ifnargin(f)>1
F(k‐1,:)=feval(f,t(k),y(k,:));
else
F(k‐1,:)=feval(f,t(k));

end


Đểgiảiphươngtrìnhtadùngchươngtrình
cthamming.m:

clearall,clc
a=0;
b=1;
y=@f1;
ya=[011]ʹ;
n=10;
tic
[t,y]=hamming(y,a,b,ya,n);
toc
plot(t,y)


§9.PHƯƠNGPHÁPMILNE
 Quátrìnhtínhtoánnghiệmđượcthựchiệnquababước:
‐Tínhgầnđúngy
m+1theocôngthức(dựđoán):

()
+− −−
′′ ′
=+ −+
(1)
m1 m3 m2 m1 m
4h

3
(3)
Taxâydựnghàm
milne()đểthựchiệnthuậttoántrên:

function[t,y]=milne(f,to,tf,y0,n)
h=(tf‐to)/n;
y(1,:)=y0ʹ;
ts=to+3*h;
[t,y]=rungekutta(f,to,ts,y0,3);
t=[t(1:3)ʹt(4):h:tf]ʹ;
fori=2:4

384
ifnargin(f)>1
F(i‐1,:)=feval(f,t(i),y(i,:));
else
F(i‐1,:)=feval(f,t(i));
end
end
fori=4:n
p=y(i‐3,:)+(4*h/3)*(2*F(1,:)‐F(2,:)+2*F(3,:));%Pt.(1)
ifnargin(f)>1
F(4,:)=f(t(i+1),
p);%Pt.(2)
else
F(4,:)=f(t(i+1));
end
y(i+1,:)=y(i‐1,:)+(h/3)*(F(2,:)+4*F(3,:)+F(4,:));%Pt.(3)
F(1,:)=F(2,:);