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

Phan Quang Tuyển

MỘT SỐ THUẬT TOÁN RUNGE - KUTTA VỚI
BƯỚC LƯỚI THAY ĐỔI GIẢI
MỘT LỚP PHƯƠNG TRÌNH VI PHÂN ĐẠI SỐ

LUẬN VĂN THẠC SĨ KHOA HỌC

Hà Nội - 2019


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

Phan Quang Tuyển

MỘT SỐ THUẬT TOÁN RUNGE - KUTTA VỚI
BƯỚC LƯỚI THAY ĐỔI GIẢI
MỘT LỚP PHƯƠNG TRÌNH VI PHÂN ĐẠI SỐ

Chuyên ngành: Toán ứng dụng
Mã số: 8460112.01
LUẬN VĂN THẠC SĨ KHOA HỌC
NGƯỜI HƯỚNG DẪN KHOA HỌC
PGS. TSKH. Vũ Hoàng Linh

Hà Nội - 2019



13

1.1.1

Khái niệm và phân loại phương trình vi phân đại số . . . .

13

1.1.2

Chỉ số của phương trình vi phân đại số . . . . . . . . . . .

15

Phương pháp Runge-Kutta cho phương trình vi phân thường . .

18

1.2.1

Phương pháp Runge-Kutta tổng quát . . . . . . . . . . . .

19

1.2.2

Sự hội tụ và tính ổn định của phương pháp Runge-Kutta .

20


25

2.2

Trường hợp phương trình vi phân đại số không có tính lạ . . . . .

27

2.2.1

Phân tích bài toán . . . . . . . . . . . . . . . . . . . . . . . .

28

2.2.2

Phương pháp Runge-Kutta nửa hiện . . . . . . . . . . . . .

30

2.3

Trường hợp phương trình vi phân đại số không có tính lạ và có
cấu trúc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.3.1


3

44

Phương pháp Runge-Kutta với bước lưới thay đổi giải phương trình vi
phân đại số

48

3.1

Phương pháp nhúng . . . . . . . . . . . . . . . . . . . . . . . . . .

48

3.1.1

Phương trình vi phân đại số không có tính lạ . . . . . . . .

50

3.1.2

Phương trình vi phân đại số không có tính lạ và có cấu trúc

51

3.2

Thử nghiệm số . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Rm , (Cm )

Không gian véc tơ m chiều trên R, (C)

Rm1 ,m2

Không gian ma trận thực cỡ m1 × m2

C p (I, Rm )

Không gian các hàm véc tơ m chiều khả vi liên tục cấp p

C p (I, Rm1 ,m2 )

Không gian các hàm ma trận cỡ m1 × m2 khả vi liên tục cấp p

Ik

Ma trận đơn vị cấp k

rank( A)

Hạng của ma trận A

Const

Hằng số nào đó

O(hk )


Phương pháp Runge-Kutta nửa hiện (Half explicit Runge-Kutta)

IRK

Phương pháp Runge-Kutta ẩn (Implicit Runge-Kutta)

HEOL

Phương pháp một chân nửa hiện (Half explicit One - leg )

HELM

Phương pháp đa bước nửa hiện (Half explicit linear multistep)

7


LỜI MỞ ĐẦU
Trong thực tế, chúng ta gặp rất nhiều bài toán trong các lĩnh vực khoa học kĩ
thuật như cơ học, hóa học, hệ mạch điện, lý thuyết điều khiển, động lực học chất
lỏng, v.v. được mô hình hóa dưới dạng một hệ hỗn hợp các phương trình vi phân
kết hợp với các ràng buộc đại số. Các hệ đó được gọi là phương trình vi phân
đại số (PTVPĐS, DAEs). PTVPĐS có dạng tổng quát
F (t, x, x ) = 0,

(0.0.1)

trong đó t ∈ I = [0, T ], F : I × Rm × Rm → Rn , m, n ∈ N. Nếu ma trận Jacobi
của F theo x không suy biến, theo định lý hàm ẩn, từ phương trình (0.0.1) ta có


(0.0.2)

x2 + y2 − l 2 = 0.
Bằng cách đặt biến mới u = x ,

v = y , phương trình (0.0.2) trở thành
x − u = 0,
y − v = 0,
mu + 2λx = 0,

(0.0.3)

mv + mg + 2λy = 0,
x2 + y2 − l 2 = 0.
Đây chính là một PTVPĐS có chỉ số 3 với các biến là x, y, u, v, λ.
Giả sử, ta đi giải bài toán giá trị ban đầu (BTGTBĐ) (0.0.1) với điều kiện đầu
x (0) = x0 , x0 ∈ Rm . Sự tồn tại và duy nhất nghiệm của BTGTBĐ (0.0.1) phụ
thuộc vào điều kiện ban đầu x0 .
Trong ví dụ (0.0.1) để hệ (0.0.3) có nghiệm thì điều kiện ít nhất ta cần có là
x02 + y20 = l02 . Không những vậy, điều kiện ban đầu của PTVPĐS còn có thể liên
quan đến đạo hàm của các ràng buộc tại thời điểm ban đầu, xem [3].
Các PTVPĐS xuất hiện từ các bài toán thực tế thường là các hệ rất phức tạp,
không có hy vọng giải đúng, trong khi nhiều trường hợp chúng ta chỉ cần biết
thông tin về nghiệm số hoặc nghiệm gần đúng với mức độ chính xác nhất định
nào đó.
Việc nghiên cứu giải số cho PTVPT đã phát triển từ lâu, trong khi việc nghiên
cứu lý thuyết và các phương pháp số của PTVPĐS mới phát triển mạnh trong
khoảng hơn 30 năm trở lại đây. Các phương pháp số cho PTVPĐS đều được mở
9

Tác giả V.H. Linh và V. Mehrmann [7] đã nghiên cứu tính chất của bài toán

(0.0.4) và đề xuất các phương pháp một chân nửa hiện (HEOL), phương pháp
đa bước nửa hiện (HELM), phương pháp Runge-Kutta nửa hiện (HERK) để giải
hiệu quả các PTVPĐS (0.0.4) không có tính cương.
Khi nghiên cứu một lớp PTVPĐS có cấu trúc dạng
f (t, x, E(t) x ) = 0,

(0.0.6)

g(t, x ) = 0,
trên đoạn [0, T ], trong đó E ∈ C1 (I, Rm1 ,m ) và các hàm f = f (t, u, v) : I × Rm ×
Rm1 → Rm1 , g = g(t, u) : I × Rm → Rm2 , (m, m1 , m2 ∈ N, m = m1 + m2 ) đủ trơn
10


và có các đạo hàm riêng bị chặn. Giả sử BTGTBĐ có nghiệm duy nhất x (t) và
fv E
gu

không suy biến dọc theo nghiệm x(t).

(0.0.7)

Đây là một lớp các PTVPĐS nằm trong dạng PTVPĐS không có tính lạ (0.0.4).
Hai tác giả V.H. Linh và N.D. Trường [8] đã nghiên cứu và đưa ra các phương
pháp Runge-Kutta nửa hiện (HERK) và phương pháp Runge-Kutta ẩn (IRK) để
tìm nghiệm số và đánh giá sai số, cấp chính xác.
Phương pháp Runge-Kutta nửa hiện đã được kế thừa từ phương pháp RungeKutta cho PTVPT được áp dụng cho lớp bài toán (0.0.4) và (0.0.6). Phương pháp
HERK được các tác giả V.H. Linh và V. Mehrmann [7], V.H. Linh và N.D. Trường

Chương 1
Giới thiệu
Trong chương này, chúng ta sẽ nhắc lại một số khái niệm cơ bản và các kiến
thức bổ trợ sẽ được sử dụng trong luận văn. Phần đầu tiên, chúng ta sẽ giới
thiệu về PTVPĐS. Phần thứ hai, tôi sẽ trình bày tóm tắt về phương pháp RK cho
PTVPT. Phần cuối cùng của chương 1, tôi sẽ trình bày việc đánh giá sai số và lựa
chọn bước đi bằng phương pháp nhúng áp dụng cho PTVPT.

1.1. Phương trình vi phân đại số
1.1.1. Khái niệm và phân loại phương trình vi phân đại số

Phương trình vi phân đại số dạng tổng quát là phương trình
(1.1.1)

F (t, x, x ) = 0,

trong đó t ∈ I = [0, T ], F : I × Rm × Rm → Rn , m, n ∈ N, nếu ma trận Jacobi

∂F
∂x

suy biến.
Ví dụ 1.1.1. Xét hệ
x1 − x1 + 1 = 0,

(1.1.2)

x1 x2 + 2 = 0,
Ta có
F (t, x, x ) =

x1
x2 . Do đó, hệ (1.1.2) là một PTVPĐS.
Nhận xét 1.1.1. Trong ví dụ trên ta thấy đạo hàm x2 không xuất hiện. Giải phương
trình đầu tiên của (1.1.2), ta được x1 = x1 + 1. Thay x1 vào phương trình thứ hai, thì
phương trình (1.1.2) sẽ được viết lại là
x1 = x1 + 1,

(1.1.3)

( x1 + 1) x2 + 2 = 0.
Ta thấy phương trình đầu tiên của (1.1.3) là một phương trình vi phân, trong khi đó
phương trình thứ hai là một phương trình đại số. Như vậy, nói một cách dễ hiểu, thì một
PTVPĐS sẽ bao gồm các phương trình vi phân kết hợp với các ràng buộc đại số.
Để thảo luận về câu hỏi sự tồn tại và duy nhất nghiệm của PTVPĐS, chúng ta
xét BTGTBĐ (1.1.1) với điều kiện x (t0 ) = x0 , t0 ∈ I, x0 ∈ Rm .
Định nghĩa 1.1.1. ([5]).
Cho C k (I, Cn ) là một không gian véc tơ tất cả các hàm khả vi, liên tục k lần từ một
đoạn I vào không gian véc tơ phức Cn .
1. Một hàm x ∈ C1 (I, Cn ) được gọi là một nghiệm của (1.1.1) nếu nó thỏa mãn
(1.1.1) tại từng điểm.
2. Một nghiệm x của PTVPĐS(1.1.1) thỏa mãn điều kiện ban đầu được gọi là nghiệm
của BTGTBĐ.
3. Một điều kiện ban đầu x (t0 ) = x0 được gọi là tương thích với PTVPĐS (1.1.1)
nếu bài toán (1.1.1) kết hợp với điều kiện ban đầu có ít nhất một nghiệm. Khi đó,
bài toán (1.1.1) với điều kiện x (t0 ) = x0 gọi là giải được.
14


Thông thường các PTVPĐS có cấu trúc toán học tùy thuộc vào phạm vi ứng
dụng nhất định. Do đó, chúng ta có các hệ PTVPĐS phi tuyến, PTVPĐS tuyến

qua các phép lấy đạo hàm. Thước đo này dường như không phản ánh được chính
xác bản chất của PTVPĐS bởi trong đó chúng ta hầu như chỉ quan tâm tới tính
chất vi phân mà không để ý đến đặc trưng của các ràng buộc đại số. Thực tế, các
15


ràng buộc đại số đôi khi làm cho bài toán trở nên phức tạp hoặc đôi khi làm cho
bài toán trở nên đơn giản. Tiếp theo, chúng ta đề cập tới khái niệm chỉ số lạ đã
được P. Kunkel và V. Mehrmann đưa ra phản ánh được cả bản chất vi phân và
các đặc trưng của phần đại số của PTVPĐS. Để định nghĩa chỉ số lạ, chúng ta xét
hệ sau
F (t, x, x , . . . , x (

+1)

(1.1.4)

) = 0,

trong đó F có dạng

F (t, x, x , . . . , x (

+1)



) = 





và định nghĩa các ma trận Jacobi
M (t, x, x , . . . , x (

+1)

) = F ;x ,...,x( +1) (t, x, x , . . . , x (

N (t, x, x , . . . , x (

+1)

+1)

) = − F ;x (t, x, x , . . . , x (

+1)

),

), 0, . . . , 0

(1.1.5)

Giả thiết 1.1.1. ([5]). Tồn tại các số nguyên µ, a, d sao cho tập nghiệm
Lµ = {(t, x, x , . . . , x (µ+1) ) ∈ R(µ+2)n+1 | Fµ (t, x, x , . . . , x (µ+1) ) = 0}
( µ +1)

khác rỗng và mỗi (t0 , x0 , x0 , . . . , x0

các biến. Từ đó ta có thể thu được PTVPT bằng việc giải biến đại số từ các ràng
buộc và thế vào các phương trình còn lại.
Lý thuyết về chỉ số lạ cho PTVPĐS phi tuyến tổng quát (1.1.1) đã được nghiên
cứu. Trong bài báo ([5]) có thuật toán biến đổi dạng (1.1.1) về một PTVPĐS dạng
không có tính lạ (0.0.4). Chúng ta quy ước rằng khi nhắc đến chỉ số của PTVPĐS
mà không nói gì thêm thì đó là chỉ số vi phân của bài toán. Dựa vào đó, chúng
ta giới thiệu lớp PTVPĐS thường gặp có dạng:
1. PTVPĐS dạng nửa hiện chỉ số 1 (Hessenberg chỉ số 1), xem ([3])
x = f (t, x, z)

(1.1.7)

0 = g(t, x, z).
Trong đó, ma trận hàm Jacobi gz được giả thiết là không suy biến với mọi
t. Từ phương trình thứ hai của (1.1.7), theo định lý hàm ẩn ta có thể giải
ra z = φ(t, x ), thế vào phương trình thứ nhất ta thu được phương trình vi
phân đối với x là x = f (t, x, φ(t, x )).
17


Như vậy, chúng ta có thể thấy PTVPĐS (1.1.7) có chỉ số vi phân bằng 1
nhưng có chỉ số lạ bằng 0 hay dạng không có tính lạ.
2. PTVPĐS dạng nửa hiện chỉ số 2 (Hessenberg chỉ số 2), xem ([3])
x = f (t, x, z)

(1.1.8)

0 = g(t, x ).
Giả thiết rằng, ma trận Jacobi gx f z không suy biến với mọi t. Biến đại số z
không xuất hiện trong phương trình thứ 2 của (1.1.8). Từ phương trình thứ


∑ aij f (tn−1 + c j hn , Yj )

(1.2.1)

∑ bi f (tn−1 + ci hn , Yi ),

(1.2.2)

j =1
s

y n = y n −1 + h n

i =1

trong đó, hn = tn+1 − tn , Yi ≈ y(tn−1 + ci hn ) là nghiệm xấp xỉ tại điểm nấc
Ti = tn−1 + ci hn , i = 1, 2, . . . , s. Các hệ số của phương pháp RK thường được cho
dưới bảng Butcher
c

A
,
bT

(1.2.3)

với A = [ aij ]s×s , b = (b1 , b2 , . . . bs ) T , c = (c1 , c2 , . . . cs ) T , chúng ta sẽ luôn chọn
s


0

1

0
0

0 0 0
0 0 0
1
2 0 0
0 1 0

1
6

1
3

1
2
1
2

• Công thức RK 4 nấc cổ điển:

1
2

1



1.2.2. Sự hội tụ và tính ổn định của phương pháp Runge-Kutta

Phương pháp RK có thể được viết lại theo phương pháp một bước dưới dạng
y n = y n −1 + h n Ψ ( t n −1 , y n −1 , h n ),

(1.2.7)

trong đó Ψ thỏa mãn điều kiện Lipschitz theo y, từ đó suy ra phương pháp RK
0- ổn định, xem ([3]).
Việc xác định cấp chính xác cho các phương pháp RK s nấc với s > 2 không
đơn giản. Chúng ta có một kết quả về cấp chính xác và sự hội của phương pháp
RK.
Định lý 1.2.1 ([5], Theorem 5.9). Nếu các hệ số aij , bi , ci của phương pháp RK thỏa
mãn các điều kiện:
s

B( p) :

1

∑ bi cik−1 = k ,

i =1
s

C (q) :

1

trong đó 1 = (1, 1, . . . 1) T . Như vậy, hàm ổn định của phương pháp RK tổng
quát trên là R(z) = 1 + zb T ( I − zA)−1 1 và miền ÔĐTĐ của phương pháp RK là
S = {z = hn λ ∈ C : | R(z)| ≤ 1}.
20


1.3. Đánh giá sai số và lựa chọn bước đi bằng phương pháp nhúng

Một số PTVPT có thể có các nghiệm mà nó thay đổi một cách nhanh chóng
trong một khoảng thời gian và lại thay đổi một cách chậm chạp trong một khoảng
thời gian khác. Do đó, chúng ta không thể sử dụng bước đi h trong phương pháp
RK như một hằng số, thay vào đó chúng ta nên sử dụng h nhỏ khi mà nghiệm
thay đổi nhanh chóng theo thời gian. Một phương pháp số mà tự động lựa chọn
bước lưới trong mỗi bước là một phương pháp thích nghi. Một lớp các phương
pháp thích nghi cho việc giải PTVPT là các phương pháp nhúng RK.
1.3.1. Ý tưởng của phương pháp nhúng RK

Trong phần này, chúng ta sẽ tìm hiểu một số cách để ước lượng sai số và điều
chỉnh bước lưới h = hn = tn − tn−1 . Cho mỗi thành phần j của y, (1 ≤ j ≤ m), ta
sử dụng sai số tương đối (RTOL) hoặc sai số tuyệt đối (ATOL j ). Chúng ta muốn
chọn h cho mỗi j, (1 ≤ j ≤ m),

|(l j )n | ≤ f rac[ ATOL j + |(y j )n | RTOL],
ở đây frac là hệ số an toàn (frac=0.8 hoặc 0.9). Để cho thuận tiện và tránh rắc rối
cho các phương pháp RK người ta sử dụng sai số địa phương l j , 1 ≤ j ≤ m.
Ý tưởng của phương pháp nhúng là chúng ta đi tính hai nghiệm xấp xỉ yn và
yˆ n tại tn , sao cho yˆ n − yn là ước lượng sai số địa phương của hai nghiệm xấp xỉ.
Chúng ta kiểm tra bất đẳng thức |yˆ n − yn | ≤ TOL. Nếu bất đẳng thức này không
thỏa mãn thì bước lưới h sẽ bị loại và bước lưới khác h˜ sẽ được lựa chọn thay thế.
Nếu phương pháp cho việc tìm yn có cấp chính xác p thì ln (h˜ ) ≈ ch˜ p+1 . Vì vậy,


• Ví dụ đơn giản nhất là phương pháp Euler hiện được nhúng trong phương

pháp hình thang:
0 0
1 1
1

0
0
0

1
2

1
2

• Phương pháp nhúng Fehlberg 4(5) có 6 nấc sử dụng phương pháp RK cấp

4 được nhúng trong phương pháp RK cấp 5. Phương pháp nhúng Fehlberg
4(5) có hệ số sai số của nghiệm yn có cấp chính xác 4 là nhỏ nhất.
• Phương pháp nhúng Dormand-Prince có 7 nấc sử dụng phương pháp RK

cấp 4 được nhúng trong phương pháp RK cấp 5. Phương pháp nhúng
Dormand-Prince 4(5) có hệ số sai số của nghiệm yˆn có cấp chính xác 5 là
nhỏ nhất. Trong khi thực hiện bước tính tiếp theo ta thường chọn nghiệm
có giá trị chính xác cao hơn thay thế cho yn−1 . Hơn nữa, ta thấy bộ hệ số
của b T trùng với hệ số tại nấc thứ 7, do đó phương pháp nhúng DormandPrince thực chất ta thực hiện có giá như phương pháp RK có 6 nấc. Vì vậy,
phương pháp nhúng Dormand-Prince sẽ có độ chính xác tốt hơn phương

1
p +1

.

Tiếp tục quá trình giống như bước tính đầu tiên cho tới khi dừng lại. Khi đó, chúng
ta có thể in ra kết quả nghiệm chính xác, sai số của hai nghiệm xấp xỉ, sai số thực tế của
nghiệm số và nghiệm chính xác, vẽ đồ thị nghiệm chính xác so với đồ thị nghiệm số, v.v.

24


Chương 2
Phương pháp Runge-Kutta nửa hiện
giải phương trình vi phân đại số
Trong chương này, tôi sẽ trình bày chi tiết các phương pháp Runge-Kutta
nửa hiện cho 3 trường hợp: PTVPĐS dạng nửa hiện, PTVPĐS không có tính
lạ, PTVPĐS không có tính lạ và có cấu trúc. Các kết quả chính của chương này
dựa vào [7], [8].

2.1. Trường hợp phương trình vi phân đại số dạng nửa hiện chỉ số 1

Các phương pháp số cho PTVPĐS dạng nửa hiện chỉ số 1 đã được nghiên cứu
và phát triển trong những năm cuối của thể kỉ 20. Hướng nghiên cứu hết sức tự
nhiên bằng cách mở rộng các phương pháp số của PTVPT cho PTVPĐS. Trong
phần này, tôi sẽ trình bày việc rời rạc hóa bài toán PTVPĐS dạng nửa hiện chỉ số
1 bằng phương pháp RK.
Các PTVPĐS dạng nửa hiện chỉ số 1 là PTVPĐS đơn giản nhất có dạng
y (t) = F (t, y(t), z(t)),