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 6
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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Phương pháp số giải phương trình vi phân thường . . . . . . . . . . . 12
1.2.1 Các khái niệm cơ bản . . . . . . . . . . . . . . . . . . . . . . . . 12
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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3 Nghiệm số của phương trình vi phân có chậm: Phương pháp cho
phương trình vi phân thường liệu có đủ hay không? . . . . . . . . . . . 17
1.3.1 Sự thất bại về cấp chính xác của phương pháp . . . . . . . . . . 18
1.3.2 Sự thất bại về tính ổn định . . . . . . . . . . . . . . . . . . . . 20
1.3.3 Một phương pháp tốt cho các PTVPCC . . . . . . . . . . . . . 25
2 Sự tồn tại và tính chính qui của nghiệm của phương trình vi
phân có chậm 27
2.1 Vị trí của các điểm gián đoạn và sự trơn dần của nghiệm . . . . . . . 27
2.1.1 Các điểm gián đoạn gốc và thứ cấp . . . . . . . . . . . . . . . . 28
2.1.2 Chậm triệt tiêu và không triệt tiêu . . . . . . . . . . . . . . . . 31
2.1.3 Chậm bị chặn và không bị chặn . . . . . . . . . . . . . . . . . . 32

1
), , y(t − τ
n
)), t ≥ t
0
,
với τ
i
= τ
i
(t, y(t)) ≥ 0 ∀t ≥ t
0
, 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
khác biệt quan trọng nhất về tính chất định tính cũng như về khía cạnh giải số
sẽ được đề cập đến trong phần đầu của chương. Phần thứ hai của chương giới
thiệu một cách vắn tắt các khái niệm cơ bản về phương pháp số giải PTVPT,
từ đó phần thứ ba sẽ chỉ ra lí thuyết cho PTVPT là không đủ khi áp dụng cho

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 ≥ t
0
,
y(t
0
) = y
0
,
(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
khi ta cần biến đổi vế phải của (1.1) để thể hiện sự phụ thuộc của biến y

vào
các giá trị trong quá khứ của biến trạng thái y. Dạng tổng quát nhất của các

,
y
t
0
= y(t
0
+ θ) = Φ(θ),
(1.2)
trong đó Φ(θ) ∈ C biểu diễn điểm khởi tạo hoặc dữ liệu khởi tạo.
7
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
Trong tài liệu này, bài toán giá trị ban đầu (1.2) sẽ được mô tả theo một
cách thức thân thiện hơn như sau

y

(t) = f(t, y(t − τ
1
), y(t − τ
n
)), t ≥ t
0
,
y(t) = φ(t), t ≤ t
0
.
(1.3)
Tùy theo độ phức tạp của hiện tượng, các chậm τ
i

≡ 0, khi đó (1.3)
có dạng như sau

y

(t) = f(t, y(t), y(t − τ)), t ≥ t
0
,
y(t) = φ(t), t ≤ t
0
.
(1.4)
Vì với t ≥ t
0
nào đó có thể xảy ra t − τ < t
0
, 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 y
0
như đối với (1.1). Nói
chung, đạo hàm bên phải y

(t
+
0
), đó là f(t
0
, φ(t
0

y

(t) = −y(t − 1), t ≥ 0,
y(t) = 1, t ≤ 0.
(1.5)
8
Nghiệm của phương trình được miêu tả trong Hình 1.1. Vì y

(0
−
) = 0 và
y

(0
+
) = −y(−1) = −1, hàm đạo hàm y

(t) có một điểm gián đoạn tại t = 0. Đạo
hàm cấp hai y

(t) được cho bởi
y

(t) = −y

(t − 1),
và do đó nó có một điểm gián đoạn tại t = 1. Đạo hàm cấp ba được cho bởi
y

(t) = −y

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 t
0
,
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(t
0
, φ(t
0
), φ(t
0
− τ), φ

(t
0
− τ))
đượ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

dao động chỉ khi hệ thống có ít nhất hai thành phần và có thể dao động hỗn
loạn chỉ khi hệ thống có ít nhất ba thành phần (Định lí Poincaré-Bendixon), thì
các nghiệm của PTVPCC có thể có tính chất dao động và thậm chí dao động
hỗn loạn trong trường hợp vô hướng.
Ví dụ 1.1.3. Xét phương trình logistic có chậm sau
y

(t) = ay(t)(1 − y(t − 1)), (1.8)
phương trình này mô hình hóa sự thay đổi của dân số, là cải tiến của mô hình
Verhulst-Pearl y

(t) = ay(t)(1 − y(t)).
Trong khi các nghiệm của phương trình Verhulst-Pearl là đơn điệu, các
nghiệm dương của (1.8) là đơn điệu với a ∈ (0, 1/e), dao động với a ∈ [1/e, π/2)
và xấp xỉ với các quỹ đạo tuần hoàn với a > π/2 (Hình 1.3 và 1.4). (Xem [6])
Cuối cùng, ta thấy rằng sự có mặt của thành phần chậm có thể thay đổi
mạnh mẽ tính chất định tính của nghiệm bằng cách tác động đến sự ổn định
của mô hình.
10
Hình 1.3: Nghiệm của (1.8) với y(t) = 0.1 khi t ≤ 0, a = 1.4 và 0.3.
Hình 1.4: Nghiệm của (1.8) với a = 1.7 trong mặt phẳng pha.
Ví dụ 1.1.4. Xét phương trình vô hướng tuyến tính

y

(t) = λy(t) + µy(t − 1), t ≥ 0,
y(t) = −t + 1, t ≤ 0,
(1.9)
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

thành phần làm mất ổn định (destabilizer).
1.2. Phương pháp số giải phương trình vi phân thường
Để có cơ sở theo dõi được điểm khác biệt giữa các phương pháp số giải
PTVPT và PTVPCC, trong phần này ta sẽ nhắc lại các khái niệm cơ bản liên
quan đến PTVPT, một số phương pháp số tiêu biểu giải PTVPT cũng như tính
chất của các phương pháp đó. (Xem [3])
1.2.1. Các khái niệm cơ bản
Xét bài toán giá trị đầu cho PTVPT tổng quát

y

= f(t, y), 0 ≤ t ≤ b,
y(0) = c,
(1.12)
trong đó b > 0, c ∈ R
n
cho trước, f ∈ C([0, b] × R
n
→ R
n
).
Định nghĩa 1.2.1. Nghiệm y(t) được gọi là ổn định nếu với mọi ε > 0 cho
trước, tồn tại một δ > 0 sao cho với mọi nghiệm ˆy(t) khác y(t) và thoả mãn
|y(0) − ˆy(0)| ≤ δ
ta cũng có
|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 ∆ = {t

h
n
− f(t
n−1
, y
n−1
) = 0.
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
N
h
u(t
n
) ≡
u(t
n
) − u(t
n−1
)
h
n
− f(t
n−1
, u(t
n−1
))
với n = 1, . . . , N. Toán tử sai phân này thay đổi tùy theo phương pháp. Xét y
h
là một hàm lưới nhận giá trị y
n

n
)
với p là một số nguyên dương.
Đặt h = max
1≤n≤N
h
n
.
Định nghĩa 1.2.5. Phương pháp số được gọi là hội tụ cấp p nếu sai số toàn cục
e
n
, trong đó e
n
= y
n
− y(t
n
), e
0
= 0 thoả mãn
e
n
= O(h
p
)
với n = 1, 2, . . . , N.
Định nghĩa 1.2.6. Phương pháp số được gọi là ổn định – 0 nếu có các hằng số
h
0
và K sao cho với các hàm lưới x

)|

, 1 ≤ n ≤ N.
Định lý 1.2.1. Nếu phương pháp là chính xác cấp p và ổn định – 0 thì nó hội
tụ cấp p
|e
n
| ≤ K max
j
|d
j
| = O (h
p
) .
Xét phương trình thử
y

= λy
với λ ∈ C, Reλ < 0 và y
0
, y
1
, , y
N
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
|y
n
| ≤ |y

= y
n−1
+ h
n
f(t
n
, y
n
).
Cũng như phương pháp Euler hiển, người ta chứng minh được rằng phương
pháp Euler ẩn chính xác cấp 1, ổn định – 0 và do đó hội tụ cấp 1. Tuy nhiên,
không giống phương pháp Euler hiển, phương pháp Euler ẩn là ổn định – A.
(c) Phương pháp trung điểm y
n
= y
n−1
+ h
n
f(t
n−1/2
,
y
n
+y
n−1
2
).
(d) Phương pháp hình thang y
n
= y

ij
f(t
n−1
+ c
j
h, Y
j
), i = 1, 2, , s,
y
n
= y
n−1
+ h
s

i=1
b
i
f(t
n−1
+ c
i
h, Y
i
),
15
trong đó




12
. . . a
1s
c
2
a
21
a
22
. . . a
2s
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
c
s
a
s1

α
j
y
n−j
= h
k

j=0
β
j
f
n−j
,
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à
t
n−j
= t
n
− jh, y


j=1
β
j
f
n−j
trong đó
β
j
= (−1)
j−1
k−1

i=j−1

i
j − 1

γ
i
,
γ
i
= (−1)
i
1

0

−s

y
n
=hf (t
n
, y
n
) ,
trong đó
∇
0
f (t
n
, y
n
) = f (t
n
, y
n
) ,
∇
i
f (t
n
, y
n
) = ∇
i−1
f (t
n
, y


ω

n+1
(t) = f(t, ω
n+1
(t), x(t − 1)), t
n
≤ t ≤ t
n+1
,
ω
n+1
(t
n
) = y
n
,
(1.14)
trong đó
x(s) =

φ(s), s ≤ 0,
η(s), 0 ≤ s ≤ t
n
.
Công thức tích phân cho ta giá trị y
n+1
và nghiệm xấp xỉ η của (1.13) liên
tục trong [t

tính ổn định được mong đợi của phương pháp PTVPT cơ sở có thể bị phá huỷ
khi phương pháp được áp dụng cho một PTVPCC.
1.3.1. Sự thất bại về cấp chính xác của phương pháp
Để minh hoạ cho khả năng có thể mất tính chính xác, xét lớp các phương
trình tuyến tính hệ số hằng số sau

y

(t) = ay(t) −
π
2
e
a
y(t − 1), t ≥ 0,
y(t) = φ(t) = e
at
sin(
π
2
t), t ≤ 0,
(1.15)
Nghiệm của phương trình, y(t) = e
at
sin(
π
2
t), thuộc lớp C
∞
trong [−1, +∞).
Theo (1.14), với n = 0, 1, . . . , ta sẽ giải PTVPT sau


π
2
s

, s ≤ 0,
η(s), 0 ≤ s ≤ t
n
Một lớp phương pháp tốt để tích phân (1.16) được cho bởi phương pháp
trùng khớp tại ν điểm Gauss, có thể được xem là các phương pháp Runge-Kutta
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
y
n+1
= y
n
+ hf

t
n
+
h
2
,

2
− e
a
π
2
φ

t
n
+
h
2
− 1

, t
n
+
h
2
− 1 ≤ 0,
y
n
+ h

a
y
n
+y
n+1
2


t
n
+
h
2
− 1

= η

t
n
+
h
2
− (m − δ) h

=

1
2
− δ

y
n−m
+

δ +
1
2

1 +
1
2
ha

y
n
−
π
2
he
a
(
n+
1
2
)
h
sin

π
2

n +
1
2

h − 1

1 −

2
−δ
)
y
n−m
+
(
δ+
1
2
)
y
n−m+1
)
1−
1
2
ha
, 0 ≤ δ ≤
1
2
,
(
1+
1
2
ha
)
y
n

(1.19)
với n > m − δ −
1
2
.
Ta biết rằng, trong khi cấp chính xác rời rạc và đều của công thức trung
điểm (ν = 1) cho các PTVPT đều bằng 2, thì công thức trùng khớp Gauss 2-nấc
(ν = 2) lại có cấp chính xác rời rạc bằng 4 và cấp chính xác đều bằng 3.
Để kiểm tra cấp chính xác rời rạc p của các phương pháp số giải PTVP với
đầu ra liên tục thu được, xét nghiệm của (1.15) trong đoạn [0, 10] và nhớ rằng
20
nghiệm thuộc lớp C
∞
. Do đó không có điểm gián đoạn nào ảnh hưởng được đến
cấp chính xác của phương pháp. Ta kí hiệu e
h
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ố r
h
= e
h
/e
h
2
được kì vọng bằng 2
p
.
Hình 1.8: Đồ thị logarit của các tỉ số r
h

=

1 +
1
2
hλ

y
n
−
4
5
hλ
1 −
1
2
hλ
21
Hình 1.9: Nghiệm của (1.20) với λ < 0.
với n ≤ m − δ −
1
2
và
y
n+1
=



(

hλ
, 0 ≤ δ ≤
1
2
,
(
1+
1
2
hλ
)
y
n
−
4
5
hλ
((
3
2
−δ
)
y
n−m+1
+
(
δ−
1
2
)

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
y
n+1
= y
n
+
h
2
(f (t
n
, y
n
, x (t
n
− 1)) + f (t
n+1
, y
n+1
, x (t
n+1
− 1))) . (1.22)
Áp dụng (1.22) cho (1.20) với cỡ bước hằng số h = 1/ (m − δ) , m ∈ Z, m ≥
2, 0 ≤ δ < 1 và nội suy tuyến tính giữa các điểm nút cho ta
y
n+1
=
(
1+

4
5
λ(2−δ+δy
1
)
1−
1
2
hλ
với n = m − 1,
và
y
n+1
=

1 +
1
2
hλ

y
n
−
1
2
h
4
5
λ ((1 − δ) y
n−m

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.
trong đó λ(t) = −50 sin
2

2π
3

t −
1
4

, có nghiệm được vẽ trong Hình 1.12 (bên
trái), là ổn định tiệm cậ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) = y
0
, ta biết
rằng nghiệm tiến đến 0 với mọi hàm thực λ(t) ≤ 0 sao cho
t

0
λ(s)ds → −∞ khi
t → +∞. Ta cũng biết rằng phương pháp trung điểm là ổn định đại số và như
một hệ quả, nghiệm số bị chặn bởi giá trị khởi tạo |y
0
| với cỡ bước hằng số h > 0

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/2 −1/2
1 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
Y
n
1
= y
n
+
1
2
h

λ(t
n
)Y
n
1
− λ(t
n+1
)Y

n
)Y
n
1
+ λ(t
n+1
)Y
n
2
−
4
5
λ(t
n
)η(t
n
− 1) −
4
5
λ(t
n+1
)η(t
n+1
− 1)

,
y
n+1
= Y
n

n+1
− 1) = (1 − δ)y
n−m+1
+ δy
n−m+2

với n ≥ m.
Hình 1.15: Nghiệm số của phương trình (1.24) đạt được bằng phương pháp Lobatto
IIIC hai nấc với λ(t) = −50 sin
2

2π
3

t −
1
4

và m − δ = 2.