1
ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN

ĐỖ THỊ THÚY NGỌC

PHƯƠNG PHÁP SỐ
GIẢI PHƯƠNG TRÌNH VI PHÂN CÓ CHẬM

Chuyên ngành : TOÁN HỌC TÍNH TOÁN
Mã số : 60 46 30

Người hướng dẫn : PGS. TS. VŨ HOÀNG LINH

HÀ NỘI - 2012


2

Mục lục

Lời nói đầu

4

1 Giới thiệu
1.1 Một vài ví dụ so sánh phương trình vi phân có chậm và phương
trình vi phân thường . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Phương pháp số giải phương trình vi phân thường . . . . . . . .
1.2.1 Các khái niệm cơ bản . . . . . . . . . . . . . . . . . . . . .
1.2.2 Một số phương pháp số tiêu biểu giải phương trình vi phân

.

.
.
.
.

.
.
.
.
.

.
.
.
.

17
18
20
25

.
.
.
.
.

27

Phụ lục

70


4

LỜI NÓI ĐẦU
Ta đã biết phương trình vi phân thường (PTVPT) là phương trình có dạng:
y (t) = f (t, y(t)), t ≥ t0 .

Tuy nhiên, trong nhiều vấn đề, biến y còn phụ thuộc vào các giá trị trong quá
khứ của biến y . Khi đó PTVPT trên biến đổi thành phương trình sau, gọi là
phương trình vi phân có chậm (PTVPCC):
y (t) = f (t, y(t − τ1 ), ..., y(t − τn )), t ≥ t0 ,

với τi = τi (t, y(t)) ≥ 0 ∀t ≥ t0 , i = 1, ..., n, được gọi là các chậm.
Việc nghiên cứu về mặt lí thuyết cũng như phương pháp số giải PTVPT đã
thu hút được sự quan tâm của các nhà toán học trong một thời gian dài với rất
nhiều những kết quả quan trọng. Ngược lại, việc nghiên cứu đối với PTVPCC
mới được quan tâm nhiều trong thời gian từ những thập niên cuối thế kỉ 20 trở
lại đây. Sự quan tâm của các nhà toán học ứng dụng dành cho các phương pháp
số giải PTVPCC ngày càng gia tăng, thể hiện qua số lượng ngày càng nhiều các
sách chuyên khảo, bài báo và công trình nghiên cứu được công bố và đăng tải
trên các tạp chí toán học uy tín.
Mục tiêu của luận văn này là giới thiệu các phương pháp số giải PTVPCC,
đặc biệt là các phương pháp số giải PTVP với đầu ra liên tục áp dụng giải bài
toán giá trị ban đầu cho PTVP có chậm hằng số hoặc chậm chỉ phụ thuộc thời
gian. Luận văn gồm 4 chương:
Chương 1 sẽ giới thiệu PTVPCC thông qua việc so sánh với PTVPT. Những

giảng dạy và hướng dẫn nhiệt tình của các thầy các cô. Xin được cảm ơn tập
thể lãnh đạo, chuyên viên phòng Giáo dục Trung học, Sở Giáo dục và Đào tạo
Ninh Bình đã thông cảm, động viên và tạo mọi điều kiện thuận lợi về thời gian,
công việc để tác giả hoàn thành khoá học và luận văn. Cuối cùng, xin được cảm
ơn bè bạn và gia đình đã hỗ trợ, động viên và chia sẻ những khó khăn với tác
giả trong suốt thời gian học tập vừa qua.
Do thời gian và trình độ còn hạn chế, chắc chắn bản luận văn không thể
tránh khỏi những thiếu sót, tác giả rất mong nhận được sự chỉ bảo tận tình của
các thầy cô và bạn bè đồng nghiệp, tác giả xin chân thành cảm ơn!
Hà Nội, 15.11.2012
Học viên

Đỗ Thị Thuý Ngọc


6

CHƯƠNG 1
GIỚI THIỆU
Nhiều vấn đề thực tế trong vật lí, kĩ thuật, sinh học, y học, kinh tế . . . có
thể được mô hình hóa bằng một bài toán giá trị ban đầu, hoặc còn được gọi là
bài toán Cauchy, cho các PTVPT có dạng
y (t) = g(t, y(t)), t ≥ t0 ,
y(t0 ) = y0 ,

(1.1)

trong đó hàm y(t), được gọi là biến trạng thái, biểu diễn một đại lượng nào đó
tham gia vào quá trình.
Tuy nhiên, để làm cho mô hình phù hợp hơn với các hiện tượng thực tế, đôi

t, τi = τi (t) (trường hợp chậm biến thiên hoặc phụ thuộc thời gian) hoặc thậm
chí là các hàm số của t và chính y , τi = τi (t, y(t)) (trường hợp chậm phụ thuộc
trạng thái). Luận văn này sẽ tập trung sự nghiên cứu vào hai trường hợp: chậm
phụ thuộc hằng số và chậm phụ thuộc thời gian. Để đơn giản hóa về mặt kí
hiệu, hàm φ(t) được hiểu là được định nghĩa trong [ρ, t0 ], trong đó
ρ = min

1≤i≤n

min(t − τi ) .
t≥t0

Một trường hợp khá phổ biến và thú vị là khi n = 2 và τ1 ≡ 0, khi đó (1.3)
có dạng như sau
y (t) = f (t, y(t), y(t − τ )), t ≥ t0 ,
y(t) = φ(t), t ≤ t0 .

(1.4)

Vì với t ≥ t0 nào đó có thể xảy ra t − τ < t0 , sự khác nhau đầu tiên giữa các
phương trình dạng (1.1) và (1.4) đó là nghiệm của (1.4) thường phụ thuộc vào
hàm khởi tạo φ(t) hơn là phụ thuộc vào giá trị khởi tạo y0 như đối với (1.1). Nói
chung, đạo hàm bên phải y (t+
0 ), đó là f (t0 , φ(t0 ), φ(t0 − τ )), không bằng đạo hàm
−
trái y (t0 ) và do đó nghiệm y không được liên kết trơn với hàm khởi tạo φ(t) tại
điểm t0 , ở đó chỉ có tính chất C 0 -liên tục là có thể được đảm bảo. Hơn nữa, tính
không liên tục của đạo hàm sẽ lan truyền từ điểm khởi tạo t0 theo khoảng tích
phân và tạo ra các điểm gián đoạn tiếp theo mà tại đó, nghiệm càng ngày càng
trơn hơn. Như một hệ quả, thậm chí nếu các hàm f (t, y, x), τ (t, y) và φ(t) trong


trong đó hàm φ(t) được giả thiết ít nhất là C 1 -liên tục. Phương trình (1.6) được
gọi là một PTVPCC trung tính (delay differential equation of neutral type).
Như trước đã nói, hàm khởi tạo φ(t) không liên kết trơn với nghiệm y(t) tại t0 ,
nơi chỉ có duy nhất tính liên tục được đảm bảo. Điểm gián đoạn này sẽ lan
truyền thành một tập các điểm gián đoạn mà ở đó nghiệm, không giống trường
hợp không trung tính, chỉ thuộc duy nhất lớp C 0 . Do đó, trừ phi điều kiện nối
φ (t−
0 ) = f (t0 , φ(t0 ), φ(t0 − τ ), φ (t0 − τ ))

được thỏa mãn, còn không thì nghiệm của (1.6) phải được hiểu theo nghĩa tổng
quát “hầu khắp nơi”.


9

Ví dụ 1.1.2. Xét phương trình
y (t) = −y (t − 1), t ≥ 0,
y(t) = t, t ≤ 0,

(1.7)

có nghiệm được vẽ trên Hình 1.2. Vì y (0− ) = 1 và y (0+ ) = −y (−1) = −1, đạo
hàm y (t) có một điểm gián đoạn tại t = 0. Hơn nữa, vì y (t) = −y (t − 1) với mọi
t ≥ 0, đạo hàm y (t) không liên tục tại t = 1 cũng như tại t = 2, 3, . . . là các bội
số của t = 1.

Hình 1.2: Nghiệm của (1.7).

Ví dụ sau chỉ ra rằng, trong khi các nghiệm bị chặn của các PTVPT có thể


với các hệ số λ, µ là các hằng số thực. Ta biết rằng, với µ = 0, phương trình
(1.9) có dạng
y (t) = λy(t), t ≥ 0,
y(0) = 1,

(1.10)

Nghiệm của phương trình này tiệm cận tới 0 với mọi λ âm và bùng nổ với
λ dương bất kì. Hơn nữa trong trường hợp λ âm, nghiệm còn bị chặn bởi giá


11

trị khởi tạo 1. Mặt khác, với µ = 0, thành phần chậm µy(t − 1) trong (1.9) tác
động như một thành phần cưỡng bức và các tính chất đã nêu của nghiệm có thể
không còn được thoả mãn. Đặc biệt, với mọi µ > 0, tồn tại λ < 0 sao cho nghiệm
không tiệm cận tới 0 và tồn tại giá trị λ < 0 khác sao cho nghiệm tiệm cận tới
0 nhưng không bị chặn bởi giá trị khởi tạo y(0) = 1. Các tình huống này được
minh họa trong Hình 1.5 với µ = 4, λ = −3.5 và λ = −5. Cũng vậy, với λ = 0.5 và
µ = −1, thành phần chậm −y(t − 1) tác động như là một thành phần ổn định
hoá (stabilizer) của mô hình có nghiệm ổn định bất chấp tính dương của λ (xem
Hình 1.6).

Hình 1.5: Nghiệm ổn định và không ổn định của (1.9) với λ < 0.

Hình 1.6: Nghiệm ổn định của (1.9) với λ > 0.

Ví dụ 1.1.5. PTVPCC trung tính sau là một ví dụ cho thấy chậm với giá trị
nhỏ có thể tạo ra một ảnh hưởng rộng lớn

|y(t) − yˆ(t)| ≤ ε, ∀t ≥ 0.

Định nghĩa 1.2.2. Nếu có thêm điều kiện
|y(t) − yˆ(t)| → 0 khi t → ∞

thì y(t) được gọi là ổn định tiệm cận.
Chọn lưới điểm ∆ = {t0 = 0, t1 , . . . , tN = b} và đặt hn = tn − tn−1 với n =
1, . . . , N . Ta xét một phương pháp số đơn giản giải (1.12), đó là phương pháp
Euler hiển, có dạng
yn = yn−1 + hn f (tn−1 , yn−1 ).

Ta viết lại công thức trên thành
yn − yn−1
− f (tn−1 , yn−1 ) = 0.
hn


13

Cho u là một hàm bất kì định nghĩa trên lưới điểm ∆. Xét toán tử sai phân
Nh u(tn ) ≡

u(tn ) − u(tn−1 )
− f (tn−1 , u(tn−1 ))
hn

với n = 1, . . . , N . Toán tử sai phân này thay đổi tùy theo phương pháp. Xét yh
là một hàm lưới nhận giá trị yn tại mỗi điểm tn , n = 0, 1, . . . , N . Khi đó phương
pháp số cho bởi phương trình toán tử
Nh yh (tn ) = 0, n = 0, ..., N

Xét phương trình thử
y = λy

với λ ∈ C, Reλ < 0 và y0 , y1 , ..., yN là một lời giải số.


14

Định nghĩa 1.2.7. Điều kiện sau được gọi là điều kiện ổn định tuyệt đối
|yn | ≤ |yn−1 | , n = 1, 2, ..., N.

Định nghĩa 1.2.8. Miền ổn định tuyệt đối của một phương pháp số là miền
trong mặt phẳng phức z sao cho khi áp dụng phương pháp cho phương trình thử,
với z = λh nằm trong miền này, ta được một nghiệm xấp xỉ thoả mãn điều kiện
ổn định tuyệt đối.
Định nghĩa 1.2.9. Một phương pháp số được gọi là ổn định – A nều miền ổn
định tuyệt đối của nó chứa toàn bộ nửa bên trái mặt phẳng phức z = λh.

1.2.2. Một số phương pháp số tiêu biểu giải phương trình vi
phân thường
Trong phần này ta xét hai loại phương pháp tiêu biểu để giải (1.12), đó là
các phương pháp một bước và các phương pháp đa bước.
Một lớp phương pháp một bước quan trọng để giải (1.12) là các phương pháp
Runge – Kutta (R - K). Sau đây là một số phương pháp R - K đơn giản nhất.
(a) Phương pháp Euler hiển yn = yn−1 + hn f (tn−1 , yn−1 ).
Người ta chứng minh được rằng phương pháp Euler hiển chính xác cấp 1, ổn
định – 0 và do đó hội tụ cấp 1. Tuy nhiên phương pháp Euler hiển không ổn
định – A.
(b) Phương pháp Euler ẩn yn = yn−1 + hn f (tn , yn ).
Cũng như phương pháp Euler hiển, người ta chứng minh được rằng phương


s

aij , i = 1, ..., s,
j=1

s





bi = 1.
i=1

Nếu aij = 0 ∀j > i thì ta có một công thức R – K hiển. Ngược lại ta có một
công thức R – K ẩn.
Ta thể hiện phương pháp một cách thuận tiện bằng cách sử dụng bảng
Butcher như sau
c1
c2

a11 a12 . . . a1s
a21 a22 . . . a2s

cs

as1 as2 . . . ass
b1 b2 . . . b3



1/2
1

(d) Phương pháp hình thang
0
1

0
0
1/2 1/2
1/2 1/2

Ngoài các phương pháp một bước, người ta còn xây dựng các phương pháp
đa bước để giải bài toán giá trị ban đầu cho các PTVPT. Một phương pháp k bước tuyến tính tổng quát có dạng
k

k

αj yn−j = h
j=0

βj fn−j ,
j=0

trong đó các hằng số αj , βj , j = 0, . . . , k là các hệ số của phương pháp, α0 =
0; |αk | + |βk | = 0 và
tn−j = tn − jh, yn−j ≈ y(tn−j ), fn−j = f (tn−j , yn−j ) ≈ f (tn−j , y(tn−j )), j = 0, ..., k.



−s
i

γi ,

ds.

Các phương pháp Adams ẩn, còn được gọi là các phương pháp AdamsMoulton, có dạng
k

yn = yn−1 + h

βj fn−j .
j=0

Cấp của một phương pháp Adams-Moulton k -bước là p = k + 1. Các phương
pháp Adams-Moulton có các hằng số sai số nhỏ hơn, sử dụng ít hơn một bước
và có miền ổn định rộng hơn nhiều so với các phương pháp Adams-Bashforth có
cùng cấp.
Một phương pháp k -bước BDF với cấp p = k có dạng sau
k

i=1

1 i
∇ yn =hf (tn , yn ) ,
i

trong đó
∇0 f (tn , yn ) = f (tn , yn ) ,


trong đó
x(s) =

φ(s), s ≤ 0,
η(s), 0 ≤ s ≤ tn .

Công thức tích phân cho ta giá trị yn+1 và nghiệm xấp xỉ η của (1.13) liên
tục trong [tn , tn+1 ] sao cho η(tn+1 ) = yn+1 .
Một tính chất riêng của hướng tiếp cận này là trong khi phương pháp số giải
PTVPT chỉ cung cấp các giá trị xấp xỉ của nghiệm tại các điểm nút thì việc cài
đặt phương pháp số giải (1.14) có thể yêu cầu các hiểu biết về nghiệm xấp xỉ η(t)
tại một vài điểm t − 1 có thể khác các điểm nút. Do đó, nói chung, các phương
pháp số giải PTVPCC sẽ dựa trên cơ sở các mở rộng liên tục của các phương
pháp số giải PTVPT. Điều này có thể thực hiện bằng một phép nội suy hậu
nghiệm (a posteriori interpolation) của các giá trị yn được cung cấp bởi phương
pháp PTVPT rời rạc cơ sở hoặc, tốt hơn, bằng phương pháp số giải PTVP với
đầu ra liên tục, đó là các phương pháp cung cấp từng bước một xấp xỉ liên tục
của nghiệm (xem Hình 1.7). Như chúng ta sẽ thấy, sự thành công của phương
pháp số giải PTVP thường với đầu ra liên tục đạt được về mặt sự chính xác và


18

Hình 1.7: Nghiệm xấp xỉ của (1.13) đạt được bằng phương pháp số giải PTVP với đầu
ra liên tục.

tính ổn định phụ thuộc vào sự lựa chọn đặc biệt phương pháp rời rạc cũng như
là mở rộng liên tục.
Ta cũng đã chỉ ra rằng sự xuất hiện của thành phần có chậm có thể thay


19
ν -nấc cấp 2ν . Vì chúng là các phương pháp trùng khớp dựa trên sự xấp xỉ đa

thức từng khúc bậc ν , nên chúng cũng cung cấp một mở rộng liên tục η(t) có
cấp chính xác đều ν + 1. Do đó, phương pháp trùng khớp Gauss xuất hiện như
là một lớp các phương pháp liên tục hấp dẫn để tích phân các PTVPCC giống
như (1.15).
Với ν = 1, phương pháp đã được biết đến như là phương pháp trung điểm,
và với phương trình tổng quát (1.14) nó có dạng
yn+1 = yn + hf

h yn + y n+1
h
tn + ,
, x tn + − 1
2
2
2

.

(1.17)

Áp dụng (1.17) cho (1.16) với cỡ bước tích phân hằng số h = 1/(m−δ), m ≥ 2,
m nguyên, và 0 ≤ δ < 1, ta được
yn + h a yn +y2 n+1 − ea π2 φ tn + h2 − 1
yn + h a yn +y2 n+1 − ea π2 η tn + h2 − 1

yn+1 =

Tóm lại, công thức trung điểm cho (1.15) có dạng
1

yn+1 =

1 + 21 ha yn − π2 hea(n+ 2 )h sin

π
2

n+

1
2

h−1

1 − 21 ha

với n ≤ m − δ − 12 và

1
π a
1
1
 (1+ 2 ha)yn − 2 e h(( 2 −δ1)yn−m +(δ+ 2 )yn−m+1 ) , 0 ≤ δ ≤ 1 ,
2
1− ha
yn+1 =


cấp chính xác của phương pháp. Ta kí hiệu eh là sai số tuyệt đối lớn nhất tại
các điểm nút với cỡ bước tích phân h = 1/(m − δ). Bằng cách chia đôi cỡ bước,
giá trị tiệm cận của tỉ số rh = eh /e h được kì vọng bằng 2p .
2

Hình 1.8: Đồ thị logarit của các tỉ số rh , h = 1/(m − δ), m − δ = 2i ∗ k, như là một hàm
số của i, cho các nghiệm số của phương trình (1.15) với a = 1 đạt được bằng phương
pháp trùng khớp tại ν = 1 và ν = 2 các điểm Gauss. Các giá trị của m − δ được xác
định bởi k = 1 (đường nét liền) và k = 5/3 (đường nét đứt).

Với các giá trị nguyên (δ = 0) và không nguyên (δ > 0) của m − δ , tương ứng
với h là hoặc không là “ước” của chậm τ = 1, công thức trung điểm duy trì được
cấp chính xác rời rạc được kì vọng là 2. Ngược lại, công thức trùng khớp Gauss
2-nấc có cấp chính xác rời rạc là 4 hoặc 3 tuỳ theo m − δ là nguyên hay không
nguyên, trong khi cấp chính xác đều là bằng 3.

1.3.2. Sự thất bại về tính ổn định
Cũng giống như việc phân tích tính chính xác của các phương pháp số giải
PTVP với đầu ra liên tục, ta xét lớp các phương trình thử tuyến tính hệ số hằng
số sau
y (t) = λy(t) − 45 λy(t − 1), t ≥ 0,
y(t) = 1, t ≤ 0.

(1.20)

Phương trình này có nghiệm, được vẽ trên Hình 1.9 với một số giá trị của λ,
là ổn định tiệm cận với mọi λ < 0.
Công thức trung điểm (1.17), được mở rộng bởi phép nội suy tuyến tính, áp
dụng cho phương trình (1.20) có dạng
yn+1 =

1
2

(1.21)

< δ < 1,

với n > m − δ − 21 .
Ta biết rằng phương pháp trung điểm là ổn định – A, tức là với mọi phương
trình y = λy, Re(λ) < 0, phương pháp cung cấp các nghiệm số tiệm cận tới 0 với
cỡ bước tích phân h > 0 bất kì. Vì nghiệm của (1.20) ổn định tiệm cận với mọi
λ < 0, ta kì vọng tính chất đó cũng đúng cho nghiệm số đạt được bằng phương
pháp trung điểm mà không phụ thuộc vào cỡ bước h. Tuy nhiên, trong khi (1.21)
cho ta các nghiệm số ổn định của (1.20) khi m − δ là số nguyên (δ = 0), thì các
giá trị không nguyên của m − δ có thể tạo ra các nghiệm số không ổn định. Hình
1.10 thể hiện các nghiệm số của (1.20) với λ = −50, đạt được bằng cách áp dụng
(1.21) với m − δ = 10 và m − δ = 12,5. Mặc dù trong trường hợp sau ta sử dụng
một cỡ bước tích phân nhỏ hơn nhưng nghiệm đạt được lại không ổn định.
Bây giờ ta xét một phương pháp ổn định – A rất phổ biến khác, gọi là phương


22

Hình 1.10: Nghiệm số của phương trình(1.20) với λ = −50 đạt được bằng phương pháp
trung điểm với h = 1/(m − δ), m − δ nguyên hoặc không nguyên.

pháp hình thang, áp dụng cho phương trình tổng quát (1.14), như sau
yn+1 = yn +

h

cấp các nghiệm ổn định với các cỡ bước tích phân không phải là “ước” của chậm.
Hình 1.11 minh hoạ tính chất nghiệm của (1.23) với λ = −50 và m − δ = 12.5,
trong khi với các điều kiện đó, phương pháp trung điểm gặp thất bại.
Mặc dù phương pháp hình thang thể hiện tính ổn định tốt hơn phương pháp
trung điểm khi áp dụng cho các PTVPCC tuyến tính hệ số hằng số, chúng có
thể là không đủ khi áp dụng cho lớp rộng hơn các PTVPCC tuyến tính hệ số
biến thiên, thậm chí với các giá trị m − δ nguyên. Để minh hoạ cho sự khác
nhau này giữa các phương pháp khi áp dụng cho các PTVPT và PTVPCC, xét
phương trình
y (t) = λ(t)y(t) − 45 λ(t)y(t − 1), t ≥ 0,
y(t) = t + 1, t ≤ 0,

(1.24)


23

Hình 1.11: Nghiệm số của phương trình (1.20) với λ = −50 đạt được bằng phương
pháp hình thang với m − δ = 12.5.
1
trong đó λ(t) = −50 sin2 2π
3 t− 4
trái), là ổn định tiệm cận.

, có nghiệm được vẽ trong Hình 1.12 (bên

Hình 1.12: Nghiệm của phương trình (1.24) với λ(t) ≤ 0 nào đó.

Với phương trình tuyến tính hệ số biến thiên y (t) = λ(t)y(t), y(0) = y0 , ta biết
t

với cỡ bước h = 0.505, tương ứng với m − δ = 1.98, được vẽ trên Hình 1.14 mà
trong đó, sau một vài dao động giả, chúng trở nên ổn định tiệm cận.

Hình 1.13: Nghiệm số của phương trình (1.24) với λ(t) ≤ 0 nào đó đạt được bằng
phương pháp trung điểm (bên trái) và phương pháp hình thang (bên phải) với m−δ = 2.

Hình 1.14: Nghiệm số của phương trình (1.24) với λ(t) ≤ 0 nào đó đạt được bằng
phương pháp trung điểm (bên trái) và phương pháp hình thang (bên phải) với m − δ =
1.98.


25

1.3.3. Một phương pháp tốt cho các PTVPCC
Để đạt được một xấp xỉ cấp hai cho một lớp các bài toán ổn định bao gồm
cả (1.24) là một bài toán ổn định với mọi cỡ bước, ta có thể lựa chọn phương
pháp Lobatto IIIC hai nấc với nội suy tuyến tính giữa các điểm nút. Ta có bảng
Butcher của phương pháp này
0
1

1/2 −1/2
1/2 1/2
1/2 1/2

Áp dụng phương pháp cho (1.24) với cỡ bước tích phân hằng số h = 1/ (m − δ) , m ∈
Z, m ≥ 2, 0 ≤ δ < 1, ta được
Y1n = yn + 12 h λ(tn )Y1n − λ(tn+1 )Y2n − 54 λ(tn )η(tn − 1) + 45 λ(tn+1 )η(tn+1 − 1) ,
Y2n = yn + 12 h λ(tn )Y1n + λ(tn+1 )Y2n − 45 λ(tn )η(tn − 1) − 45 λ(tn+1 )η(tn+1 − 1) ,
yn+1 = Y2n ,