Bài giảng Chuyên đề Phương pháp tính - Chương 6: Nghiệm gần đúng của hệ phương trình vi phân thường

Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
Chương 6 NGHIỆM GẦN ĐÚNG CỦA HỆ 
PHƯƠNG TRÌNH VI PHÂN THƯỜNG 
6.1 Mở đầu 
 Nhiều bài toán khoa học kỹ thuật có phương trình chỉ đạo là (hệ) 
phương trình vi phân thường cùng với điều kiện biên và điều kiện ban đầu. 
Nghiệm đúng của chúng thường chỉ áp dụng cho một số lớp bài toán rất hạn 
chế; đa số các bài toán là phải tìm nghiệm gần đúng. 
Có hai loại bài toán là: 
(i) Bài toán Cauchy hay còn gọi là bài toán giá trị ban đầu, bao gồm 
(hệ) phương trình vi phân và điền kiện ban đầu của bài toán. 
(ii) Bài toán biên, bao gồm (hệ) phương trình vi phân và điều kiện biên 
Để giải gần đúng các bài toán nầy có hai phương pháp là: 
(a) Phương pháp giải tích: tìm nghiệm gần đúng dưới dạng biểu thức 
như phương pháp xấp xỉ liên tiếp Picard, phương pháp chuổi 
nguyên, phương pháp tham số bé, 
(b) Phương pháp số: tìm nghiệm gần đúng bằng số tại các điểm rời rạc; 
nó còn chia ra phương pháp một bước (như phương pháp Euler, 
Runghe-Kutta,) và phương pháp đa bước (Adams,); Với 
phương pháp một bước tính nghiệm gần đúng yi thông qua yi-1 còn 
với phương pháp đa bước yi tính được thông qua nhiều bước trước 
đó: yi-1, yi-2, yi-3, 
6.2 Nghiệm gần đúng của bài toán Cauchy đối với phương trình vi 
phân thường 
Giả sử ta cần giải bài toán Cauchy: 

=
=
00 y)x(y
)y,x(f'y
 (6.2.1) 
Giả sử rằng trong miền ta xét, hàm f(x,y) có các đạo hàm riêng liên 
tục đến cấp n, khi đó nghiệm cần tìm sẽ có các đạo hàm riêng liên tục đến 
cấp n + 1, và do đó ta có thể viết : 
)(
)!1(
)(
......"
!2
)(
)()(
1
0
)1(
0
1
0
0
2
0
,
0000
++
+
−++
−++−
+−=−=∆
nn
n
o
xxy
n
xxyxx
yxxyxyy
θ (6.2.2) 
 Ký hiệu x - x0 = h, với h đủ bé ta có thể bỏ qua 0(|x – x0|n+1). 
1
0 )(
+− nxxθ
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 33 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
Từ (6.2.2) ta có: ∆y0 = y(x0+h) - y0 + hy’0 + 
)1n(
0
1n
0
2
y
)!1n(
h.........."y
!2
h ++
+++ (6.2.3) 
Để tính (6.2.3) ta lần lượt tính từ (6.2.1): 
 y’0 = f(x0,y0) = f0 , y”0 = y
f
f
x
f 0
0
0
∂
∂+∂
∂
 , 
Nói chung ta có: ∑
= − ∂∂
∂=


∂
∂+∂
∂ n
0K
KKm
m
KK
m
m
yx
ufCu
y
f
x 
Vậy ta tính được: y(x) ≅ ∑
=
n
0K
K
0
)K(
!K
h)x(y 
Trong thực tế cách tính nầy ít dùng vì cồng kềnh; ta sẽ xét các phương 
pháp giải khác đơn giản hơn. 
6.2.1 Phương pháp xấp xỉ liên tiếp Pica 
Một trong những phương pháp giải tích giải gần đúng phương trình vi 
phân (6.2.1) là phương pháp xấp xỉ liên tiếp Pica. 
Mục đích của phương pháp này là xây dựng nghiệm cần tìm là y= y(x) 
Từ (6.2.1) ta có: ∫∫∫ =−⇒= x
x
x
x
0
x
x 000
dt)y,t(f)x(y)x(ydt)y,t(fdy
Hay: (6.2.4) ∫+= x
x
0
0
dt)y,t(fy)x(y
Giả sử f(x,y) là hàm liên tục theo x,y và y
f
∂
∂
 < K. 
Để tìm xấp xỉ liên tiếp, trong (6.2.4) thay y bằng y0, ta có xấp xỉ thứ 
nhất: 
 , ∫+= x
x
001
0
dt)y,t(fyy
Tương tự có xấp xỉ thứ hai: ∫+= x
x
102
0
dt)y,t(fyy
Tổng quát, ta có: , với n = 1,2,3, ∫ −+= x
x
1n0n
0
dt)y,t(fyy
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 34 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
Như vậy ta sẽ có: ∫ −+=≈ x
x
1n0n
0
dt)y,t(fy)x(y)x(y
 )x(y)x(ylim
n
n =∞→ 
Sai số: 
!n.K
)KC(M)x(y)x(y
n
n ≤− , trong đó )y,x(f = M 
Với: 0xx − < a ≤ ∞, 0yy− < b ≤ ∞ , thì C = min 


M
b,a 
Ta có: 
(i) y
f
∂
∂
 > 0 và f(x,y0) > 0 thì: y0 < y1 < y2 < . . . < yn < y(x) 
(ii) y
f
∂
∂
 > 0 và f(x,y0) y1 > y2 > . . . > yn > y(x) 
 Trong hai trường hợp nầy ta có dãy xấp xỉ 1 phía. 
(iii) y
f
∂
∂
 < 0 các xấp xỉ Pica lập thành các xấp xỉ 2 phía. 
6.2.2 Phương pháp Euler 
x 
y 
O xo x1 x2 x3 
Ao A1 
A2
A3 
y=f(x)
` 
Trước hết chia đọan [xo, X] thành n đọan nhỏ: 
xi=xo+ih, với i = 0,1,2,....,n 
 n
)xX(h o−= 
Đi xây dựng công thức, dùng khai triển Taylor hàm y=f(x) tại xi ta có: 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 35 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
2
iii )(!2
)( )x-).(x(xy )y(x y(x) ii xx
cy −′′+′+= 
Với: ci = xi + θ(x - xi), 0 < θ < 1 
Thay x = xi+1 = xi + h, và y’(xi) = f(xi,y(xi)) 
Ta có: !2
)c(y.h ))y(x, h.f(x ) y(x )y(x i2iii1i
′′++=+ 
Khi bước chia h khá bé, số hạng cuối ≅ 0, khi thay y(xi) bằng ui ta được: 
 ui+1 = ui + hi.f(xi,ui) 
Biểu thức nầy cho phép tính ui+1 khi biết ui, với điều kiện ban đầu được cho 
là: uo = η 
Đánh giá sai số: 
Định lý: Gỉa sử L
y
f ≤∂
∂ và Ky'' ≤ , trong đó L, K là những hằng số, khi 
đó phương pháp Euler hội tụ và sai số là ei = ui - y(xi) có đánh giá: 
2
K,eM
)he(M)x(yue
)xx(L
0iii
0i =α=
α+≤−=
− 
6.2.3 Phương pháp Runghe - Kutta bậc 4 
 Xét phương trình vi phân: u’ = f(x , u) 



++=
++=
++=
=
)ku,hx(f.hk
)k5.0u,h5.0x(f.hk
)k5.0u,h5.0x(f.hk
)u,x(f.hk
3ii4
2ii3
1ii2
ii1
⇒ ui +1 = ui + )kk2k2k(6
1
4321 +++ 
Với sai số: )h(0)x(Yu 4ii =−
6.2.4 Phương pháp Adam 
Giả sử cần giải phương trình vi phân: 
 Y’ = f(x , y), với điều kiện ban đầu: y(x0) = y0 
Cho biến số thay đổi bởi bước h nào đó; xuất phát từ điều kiện ban 
đầu Y(x0) = Y0 bằng phương pháp nào đó (ví dụ: phương pháp Runghe-
Kutta bậc 4), ta tìm được 3 giá trị tiếp theo của hàm cần tìm y(x): Y1 = 
Y(x1) = Y(x0+h), Y2 = Y(x0+2h), Y3 = Y(x0 + 3h) . 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 36 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
Nhờ các giá trị x0 , x1 , x2 , x3 và Y0 , Y1 , Y2 , Y3 , ta tính được q0, q1, 
q2, q3. 
Trong đó: q0 = h.Y0’ = h.f(x0 , y0), q1 = h.f(x1 , y1), q2 = h.f(x2 , y2), 
q3 = h.f(x3 , y3), sau đó ta lập bảng sai phân hữu hạn của các đại lượng y và q 
 x y ∆y q ∆q ∆2q ∆3q ------- 
 xo yo qo 
 ∆yo ∆q0 
 x1 y1 q1 ∆2q0 
 ∆y1 ∆q1 ∆3q0 
 x2 y2 q2 ∆2q1 -------- 
 ∆y2 ∆q2 -------- 
 x3 y3 q3 ------- 
---------
--- 
---------
- 
---------
- 
---------
- 
---------
-- 
---------
- 
--------- ------- 
Biết các số ở đường chéo dưới, ta tìm ∆y3 theo công thức Adam như sau: 
 0
3
1
2
233 q..8
3q.
12
5q.
2
1qy ∆+∆+∆+=∆ 
Tiếp đó ta có: 
 Y4 = Y3 + ∆Y3 → q4 = h.f(x4, Y4) 
Sau đó viết đường chéo tiếp theo như sau: 
 ∆q3 = q4 - q3 , ∆2q2= ∆q3 - ∆q2 , ∆3q1 = ∆2.q2 - ∆2.q 
Đường chéo mới cho phép ta tính ∆Y4 : 
∆Y4 = q4 + 1/2∆q3 + 5/12∆2q2 + 3/8∆3q1 
Vì vậy ta có: Y5 = Y4 + ∆Y4 . . . . . 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 37