Luận văn Thạc sĩ Toán học: Tính ổn định và tính co của các phương pháp Runge

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

NGUYỄN THỊ HIÊN

TÍNH ỔN ĐỊNH VÀ TÍNH CO CỦA CÁC PHƯƠNG PHÁP RUNGE-KUTTA

Chuyên ngành: Toán ứng dụng

Mã số: 60 46 01 12

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

NGƯỜI HƯỚNG DẪN KHOA HỌC: PGS.TS. VŨ HOÀNG LINH

Hà Nội-2014

Lời cảm ơn

Trước khi trình bày nội dung chính của luận văn, tôi xin cảm ơn Ban chủ nhiệm khoa Toán - Cơ - Tin học cùng toàn thể các thầy giáo, cô giáo trong khoa Toán - Cơ - Tin học, phòng Sau đại học, trường Đại học Khoa học Tự nhiên - Đại học Quốc gia Hà Nội, đã giảng dạy tận tình và tạo điều kiện thuận lợi để tôi hoàn thành tốt luận văn.

Đặc biệt, tôi xin gửi lời cảm ơn chân thành nhất tới thầy giáo PGS.TS Vũ Hoàng Linh, người đã trực tiếp hướng dẫn và tận tình chỉ bảo tôi trong suốt quá trình tôi học tập và thực hiện luận văn.

Nhân dịp này, tôi cũng xin cảm ơn gia đình đã luôn ủng hộ và động viên

trong suốt thời gian tôi học tập.

Cuối cùng, tôi xin cảm ơn tất cả các bạn, các anh, các chị trong lớp cao học Toán khóa 2011 - 2013, đặc biệt là các anh chị chuyên ngành Toán ứng dụng khóa 2010 - 2012 và khóa 2011 - 2013 đã tận tình giúp đỡ và động viên tôi trong quá trình học tập.

Xin chân thành cảm ơn!

Hà Nội, ngày 16 tháng 01 năm 2014

Học viên

Nguyễn Thị Hiên

Mục lục

Mở đầu

Bảng ký hiệu

1 Các khái niệm cơ bản

5 5 1.1 Các phương pháp Runge-Kutta . . . . . . . . . . . . . . . . 1.2 Xây dựng các phương pháp Runge-Kutta ẩn . . . . . . . . 11 1.3 Áp dụng các phương pháp Runge-Kutta giải bài toán cương 18 . . . . . . . . . . . . . . . . . . . . . . . . 21 1.4 Các loại chuẩn

2 Tính co cho bài toán tuyến tính

26 2.1 Chuẩn Euclid (Định lý von Neumann) . . . . . . . . . . . . 29 2.2 Hàm tăng trưởng sai số với bài toán tuyến tính . . . . . . . 30 2.3 Bài toán với nhiễu phi tuyến nhỏ . . . . . . . . . . . . . . . 33 2.4 Tính co trong (cid:107).(cid:107)∞ và (cid:107).(cid:107)1 . . . . . . . . . . . . . . . . . . 37 2.5 Hệ số ngưỡng . . . . . . . . . . . . . . . . . . . . . . . . . 39

3 Tính ổn định B và tính co

42 3.1 Điều kiện Lipschitz một phía . . . . . . . . . . . . . . . . . 42 3.2 Ổn định B và ổn định đại số . . . . . . . . . . . . . . . . . 43 3.3 Một vài phương pháp Runge-Kutta ẩn ổn định đại số . . . . 46 3.4 Ổn định AN . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.5 Các phương pháp Runge-Kutta khả quy . . . . . . . . . . . 51 3.6 Định lý về sự tương đương giữa ổn định B và ổn định đại

số với các phương pháp S-bất khả quy . . . . . . . . . . . . 53 3.7 Hàm tăng trưởng sai số . . . . . . . . . . . . . . . . . . . . 56 3.8 Tính toán ϕB (x) . . . . . . . . . . . . . . . . . . . . . . . 58 Kết luận . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Tài liệu tham khảo . . . . . . . . . . . . . . . . . . . . . . . . 64

Mở đầu

• Chương 1: Các khái niệm cơ bản

Trong khoa học và kĩ thuật ta thường gặp rất nhiều bài toán liên quan tới việc giải phương trình vi phân. Có rất nhiều trường hợp nghiệm giải tích của các bài toán này là không thể tìm được. Chính vì vậy các nhà toán học đã tìm kiếm nhiều phương pháp số khác nhau để giải các bài toán trên. Trong các phương pháp số, phương pháp Runge-Kutta có nhiều tính chất ưu việt và được sử dụng rộng rãi. Luận văn trình bày về tính ổn định và tính co của các phương pháp Runge-Kutta. Xuất phát từ điều kiện ổn định tuyệt đối |yn| ≤ |yn−1| của bài toán y(cid:48) = λy, ta mở rộng đến khái niệm "tính co" khi xét bài toán tuyến tính y(cid:48) = Ay, tiếp đến là các khái niệm tính ổn định B và ổn định đại số khi xét bài toán phi tuyến. Trên cơ sở đó ta có thể lựa chọn ra phương pháp hữu hiệu và phù hợp nhất để giải các bài toán nảy sinh trong thực tế. Nội dung luận văn được tham khảo chính từ tài liệu [2] và [3]. Bố cục của luận văn bao gồm 3 chương:

• Chương 2: Tính co của bài toán tuyến tính

Luận văn trình bày các khái niệm cơ bản về phương pháp Runge-Kutta, cách xây dựng phương pháp Runge-Kutta ẩn, cùng với các kiến thức bổ trợ cho Chương 2 và Chương 3.

• Chương 3: Tính ổn định B và tính co

Luận văn trình bày các khái niệm và định lý liên quan đến tính co khi xét bài toán tuyến tính.

Luận văn trình bày khái niệm ổn định B, ổn định đại số, ổn định AN và mối quan hệ giữa các khái niệm ổn định của các phương pháp Runge-Kutta khi xét bài toán phi tuyến.

Do thời gian thực hiện luận văn không nhiều nên trong luận văn không tránh khỏi những hạn chế và sai sót. Tác giả mong nhận được sự góp ý và những ý kiến phản biện của quý thầy cô và bạn đọc.

Bảng ký hiệu

A ⊗ I

Tích tensor.

B (p), C (η), D (ζ) Bộ điều kiện về cấp chính xác.

Tập số phức.

Không gian vectơ phức n chiều.

Ma trận đơn vị.

Hàm ổn định với bài toán y(cid:48) = λ (x) y.

K (Z)

Đa thức trực giao Legendre.

Pk (x)

Xấp xỉ Padé.

Pkj (z)

Hàm ổn định của phương pháp.

R (z)

Tập số thực.

Không gian vectơ thực n chiều.

Miền ổn định.

µ (A)

Chuẩn logarit của ma trận A.

Hằng số Lipschitz một phía.

Hàm tăng trưởng sai số khi xét bài toán phi tuyến.

ϕB (x)

Hàm tăng trưởng sai số khi xét bài toán tuyến tính.

ϕR (x)

(cid:37)

Hệ số ngưỡng.





Chuyển vị của vectơ b =

bT = (b1, ..., bs)

b1 ...

    

  .   

111 = (1, ... , 1)T

Vectơ cột với các tất cả các thành phần đều bằng 1.

Chương 1

Các khái niệm cơ bản

1.1 Các phương pháp Runge-Kutta

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

Chương này trình bày các khái niệm cơ bản về các phương pháp Runge- Kutta, sự tồn tại lời giải số của phương pháp, cách xây dựng các phương pháp Runge-Kutta ẩn cùng với các kiến thức bổ trợ cho Chương 2 và Chương 3. Nội dung của chương này chỉ phát biểu các khái niệm và các kết quả phục vụ cho các chương sau. Chứng minh chi tiết của các kết quả trong chương này có thể tham khảo tại [2], [3] và [5].

Phương pháp Runge-Kutta thuộc lớp các phương pháp số một bước, được đưa ra bởi hai nhà toán học người Đức là Carl Runge (1856 - 1927) và Wilhelm Kutta (1867 - 1944).

Trước hết ta xét bài toán Cauchy của phương trình vi phân cấp một có dạng

y(cid:48) = f (t, y) , y ∈ Rn, f : R × Rn → Rn, y (t0) = y0.

(1.1)

i = 1, ..., s

Yi = yn−1 + h

aijf (tn−1 + cjh, Yj)

Định nghĩa 1.1 (xem [5]). Phương pháp Runge-Kutta s nấc cho hệ phương trình vi phân (1.1) có thể viết dưới dạng:

yn = yn−1 + h

bif (tn−1 + cih, Yi)

s (cid:80) j=1 s (cid:80) i=1

(1.2)

bi = 1. Thông thường, ta

i=1 thỏa mãn

i=1 ; {aij}s

i,j=1 ; {bi}s

s (cid:80) i=1

Trong đó Y1, ..., Ys là các giá trị nấc xấp xỉ y tại ti = tn−1 + cih (ti là các điểm nấc). Bộ hệ số: {ci}s

aij (i = 1, ..., s).

s (cid:80) j=1

chọn để có ci =

• Nếu aij = 0 với i ≤ j thì phương pháp là phương pháp Runge-Kutta hiển

• Nếu aij = 0 với i < j và có ít nhất một aii (cid:54)= 0 thì phương pháp là phương

(ERK).

• Nếu aij = 0 với i < j và aii = γ với i = 1, ..., s thì phương pháp là phương

pháp Runge-Kutta ẩn đường chéo (DIRK).

• Các trường hợp còn lại gọi là phương pháp Runge-Kutta ẩn (IRK).

pháp ẩn đường chéo đơn (SDIRK).

Để dễ dàng hình dung về phương pháp Runge-Kutta, Butcher đã đưa bộ hệ

c1 c2 ... cs

· · · · · · . . . · · · · · ·

a11 a21 ... as1 b1

a12 a22 ... as2 b2

a1s a2s ... ass bs

Bảng 1.1: Bảng Butcher.

số của phương pháp vào bảng sau:

(a) Euler hiển

0 1

0 1 2

(b) Hình thang hiển 0 1 1 2

0 0 1 2

Bảng 1.2: Một số công thức ERK.

Ví dụ 1.1. Một số công thức ERK cơ bản

(a) Euler ẩn

(b) Hình thang ẩn

1 1

1 2

1 2 1

0 1 2 1 2

Bảng 1.3: Một số công thức IRK.

Ví dụ 1.2. Một số công thức IRK cơ bản

Sự tồn tại lời giải số của phương pháp

i = 1, 2, ..., s ta thu được

ki = f (t0 + cih, y0 + h

aijkj)

s (cid:80) j=1

Xét công thức (1.2) trong trường hợp n = 1, nếu ta đặt ki = f (t0 + cih, Yi) với

y1 = y0 +

biki.

s (cid:80) i=1

(1.3)

Để xác định lời giải số y1 của phương pháp, trước hết ta cần xác định các giá trị ki từ hệ phương trình chứa ki cho bởi (1.3). Nói chung, đây là hệ phương trình phi tuyến nên trong nhiều trường hợp có thể các ki tồn tại không duy nhất. Do đó, không tồn tại lời giải số của phương pháp. Định lý sau đây cho ta điều kiện để tồn tại lời giải số của phương pháp Runge-Kutta ẩn (1.3).

h <

1 L max ti

, |aij|

Định lý 1.1 (xem [2]). Cho hàm f : R × Rn → Rn là hàm liên tục và thỏa mãn điều kiện Lipschitz theo biến y với hằng số L. Nếu

(cid:80) j

Ổn định tuyệt đối và ổn định A

thì lời giải số của phương pháp (1.3) tồn tại duy nhất với các giá trị ki được xác định duy nhất từ hệ phương trình cho bởi (1.3), các giá trị này có thể thu được bằng phương pháp lặp Newton. Hơn nữa, nếu f (t, y) là hàm khả vi, liên tục tới cấp p thì các ki (là các hàm theo biến h) cũng khả vi, liên tục cho tới cấp p.

Bằng các phương pháp khác nhau ta có thể tìm ra nghiệm số của phương trình vi phân. Tuy nhiên, nghiệm số ta tìm được liệu có tốt không, làm thế nào để có thể đánh giá được nghiệm đó. Để giải quyết vấn đề này, ta cần chỉ ra nghiệm số đó phải có tính chất tốt cho ít nhất một lớp các bài toán. Xét phương trình vi phân thử

y(cid:48) = λy, λ là hằng số, λ ∈ C, y(t0) = y0.

(1.4)

n = 1, 2, ....

Ở các phần sau, khi xét phương trình thử (1.4) ta luôn giả sử Re (λ) ≤ 0. Trong trường hợp Re (λ) ≤ 0 ta được điều kiện ổn định tuyệt đối là

|yn| ≤ |yn−1| ,

(1.5)

Giả sử y(t), (cid:101)y(t) là hai lời giải của (1.1). Từ đó ta có các định nghĩa về sự ổn định và ổn định tiệm cận với nghiệm của phương trình vi phân( 1.1).

|y (t) − (cid:101)y (t)| < ε với mọi t ≥ t0.

Định nghĩa 1.2. Nghiệm y (t) của phương trình vi phân( 1.1) gọi là ổn định nếu với mọi ε > 0, ∃δ > 0 sao cho |y (to) − (cid:101)y (to)| ≤ δ thì

|y (t) − (cid:101)y (t)| = 0.

lim t→+∞

Định nghĩa 1.3. Nghiệm y (t) của phương trình vi phân( 1.1) gọi là ổn định tiệm cận nếu nó ổn định và thỏa mãn điều kiện

Nhận xét 1.1. Đối với hệ tuyến tính, nếu nghiệm tầm thường là ổn định (ổn định tiệm cận) thì mọi nghiệm cũng ổn định (ổn định tiệm cận). Trong trường hợp này, chúng ta nói hệ ổn định (ổn định tiệm cận).

y(cid:48) = Ay ,

A ∈ Cm×m.

Xét bài toán tuyến tính

(1.6)

Định lý 1.2 (xem [5]). Nếu mọi giá trị riêng λ của ma trận A đều thỏa mãn Re (λ) ≤ 0 và các giá trị riêng có phần thực bằng 0 đều là các giá trị riêng đơn thì hệ y(cid:48) = Ay là ổn định.

Định lý 1.3 (Điều kiện cần và đủ (xem [5])). Nếu Re (λ) < 0 với mọi giá trị riêng λ của A thì hệ ổn định tiệm cận.

Để đi đến khái niệm về hàm ổn định của phương pháp Runge-Kutta, ta áp dụng phương pháp Runge-Kutta với công thức (1.2) cho phương trình thử (1.4) được lời giải số

yn = R (z) yn−1 với z = λh.

(1.7)

|R (z)| ≤ 1.

Khi đó, để điều kiện ổn định tuyệt đối (1.5) được thỏa mãn thì

(1.8)

Định nghĩa 1.4. Hàm R (z) xác định ở (1.7) được gọi là hàm ổn định của phương pháp Runge-Kutta. Tập S = {z ∈ C : |R (z)| ≤ 1} được gọi là miền ổn định tuyệt đối của phương pháp Runge-Kutta.

• Phương pháp Euler hiển có hàm ổn định R (z) = 1 + z.

• Phương pháp Euler hiển có miền ổn định tuyệt đối là hình tròn có bán kính

Ví dụ 1.3.

bằng 1, tâm −1.

s (cid:88)

i = 1, 2, ..., s,

Mệnh đề 1.1 (xem [3]). Phương pháp Runge-Kutta ẩn s nấc với

gi = y0 + h

aijf (t0 + cjh, gj),

j=1 s (cid:88)

(1.9a)

y1 = y0 + h

bjf (t0 + cjh, gj),

j=1

(1.9b)

R (z) = 1 + zbT (I − zA)−1111,

khi áp dụng cho phương trình thử y(cid:48) = λy thì có y1 = R (hλ) y0 với

111 = (1, 1, ..., 1)T .

(1.10)

i,j=1 ,

trong đó bT = (b1, ..., bs) , A = (aij)s

R (z) =

Mệnh đề 1.2 (xem [3]). Hàm ổn định của (1.9) thỏa mãn

det (cid:0)I − zA + z111bT (cid:1) det (I − zA)

STT Phương pháp

Euler ẩn

Phương pháp θ

R (z) 1 1 − z 1 + z (1 − θ) 1 − zθ

Trung điểm ẩn Hình thang ẩn

SDIRK cấp 3

(cid:27) 1 + z/2 1 − z/2 1 + z (1 − 2γ) + z2 (cid:0)1/2 − 2γ + γ2(cid:1) (1 − γz)2

Bảng 1.4: Hàm ổn định của một số phương pháp Runge-Kutta ẩn.

(1.11)

R (z) =

deg P = k, deg Q = j.

Miền ổn định tuyệt đối của các phương pháp cho trong Bảng 1.4 có thể xem trong [3]. Từ kết quả ở trên ta thấy phương pháp Runge-Kutta ẩn có hàm ổn định R (z) là hàm hữu tỉ với tử số và mẫu số có bậc ≤ s.

P (z) Q (z)

(1.12)

S ⊃ C− = {z ∈ C : Re (z) ≤ 0} ,

Định nghĩa 1.5. Một phương pháp mà miền ổn định tuyệt đối thỏa mãn

thì phương pháp đó gọi là phương pháp ổn định A (hay A-ổn định).

Phương pháp Runge-Kutta với hàm ổn định như ở (1.12) là ổn định A khi và

|R (iy)| ≤ 1 với mọi y ∈ R

chỉ khi

(1.13)

R (z) là giải tích với Re (z) < 0.

và

(1.14)

Xấp xỉ Padé của ez

Cho hàm f (z) giải tích trong lân cận của điểm 0, hai số nguyên không âm k

f (z) ≈

và j. Khi đó, ta có thể xấp xỉ hàm f (z) bởi hàm hữu tỉ

Pk (z) Qj (z)

P là đa thức bậc k, Q là đa thức bậc j và sai số của xấp xỉ tới O (cid:0)zk+j+1(cid:1). Trong trường hợp j = 0 khi đó xấp xỉ chính là khai triển Taylor của hàm f (z) tại điểm 0. Khi k = 0 thì Q (z)/P (z) là khai triển Taylor của hàm 1/f (z).

(1.15)

Tuy nhiên, với hàm bất kỳ và k, j bất kỳ thì không phải lúc nào cũng tồn tại

f (z) =

+ O (cid:0)zk+j+1(cid:1)

Pkj (z) Qkj (z)

xấp xỉ này. Khi xấp xỉ

tồn tại, thì bộ số (k, j) được gọi là xấp xỉ Padé của hàm f.

Xấp xỉ Padé của hàm mũ là trường hợp đặc biệt chúng ta cần nghiên cứu, một vài xấp xỉ Padé xấp xỉ các hàm hữu tỉ của các phương pháp quan trọng như phương pháp Gauss, phương pháp Radau, phương pháp Lobatto... Ta sẽ chỉ ra luôn tồn tại xấp xỉ Padé cho các hàm ổn định của các phương pháp này với k, j bất kỳ.

Định lý 1.4 (xem [3]). Bộ số (k, j) là xấp xỉ Padé của ez được cho bởi

Rkj (z) =

Pkj (z) Qkj (z)

(1.16)

.z +

+ ... +

Pkj = 1 +

.z +

+ ... +

Qkj = 1 −

= Pjk (−z) .

k j + k j k + j

k (k − 1) . (j + k) (j + k − 1) j (j − 1) . (k + j) (k + j − 1)

z2 2! z2 2!

k (k − 1) ...1 (j + k) ... (j + 1) (−1)jj (j − 1) ...1 . (k + j) ... (k + 1)

zk k! zj j!

trong đó

Ckj = (−1)j

k!j! (k + j)! (k + j + 1)!

Đặt

Pkj − ezQkj + Ckjzk+j+1 +

Ckjzk+j+2 = O (cid:0)zk+j+3(cid:1) .

j + 1 k + j + 2

khi đó

Sao cấp chính xác

Khái niệm sao cấp chính xác được ra đưa khi nghiên cứu về tính ổn định của

xấp xỉ Padé của hàm ez (Wanner, Hairer và Norsett 1978).

A = {z ∈ C : |R (z)| > |ez|} = {z ∈ C : |q (z)| > 1} ,

Định nghĩa 1.6. Cho tập hợp

(1.17)

R (z) ez

ở đây q (z) = . Tập A khi đó được gọi là sao cấp chính xác của R (z).

Sao cấp chính xác không so sánh |R (z)| với 1 như khi xét tính ổn định, mà nó so sánh |R (z)| với nghiệm chính xác của phương trình thử |ez| = ex (z = x + iy) và hi vọng sẽ nhận được nhiều thông tin hơn. Chúng ta luôn giả sử rằng các hệ số của hàm R (z) là số thực và sao cấp chính xác đối xứng qua trục thực. Hơn nữa, với z = iy ta có |ez| = 1, khi đó tập A là phần bù của miền ổn định tuyệt đối S trên trục ảo. Hình 1.1 và Hình 1.2 thể hiện các tập sao cấp chính xác của xấp xỉ Padé.

ez − R (z) = Czp+1 + O (cid:0)zp+2(cid:1)

Bổ đề 1.1 (xem [3]). Nếu R (z) là một xấp xỉ cấp p của ez, tức là

với C (cid:54)= 0 thì khi z → 0, A có dáng điệu như một hình ngôi sao với p + 1 cánh

π p + 1

1.2 Xây dựng các phương pháp Runge-Kutta ẩn

. sao có các góc quét bằng nhau và bằng

B (p) :

q = 1, ..., p;

bicq−1

i =

C (η) :

i = 1, ..., s, q = 1, ..., η;

aijcq−1

Phần này sẽ đề cập đến lớp các phương pháp Runge-Kutta ẩn sở hữu tính ổn định khá tốt. Việc xây dựng các phương pháp như vậy chủ yếu dựa vào các bộ điều kiện 

j =

1 q cq i q

(1.18)

j = 1, ..., s, q = 1, ..., ζ.

D (ζ) :

aij =

bicq−1 i

bj q

s (cid:80) i=1 s (cid:80) j=1 s (cid:80) i=1

(cid:1) (cid:0)1 − cq j  

Điều kiện B (p) có nghĩa là công thức cầu phương (bi, ci) có cấp p. Sự quan trọng của hai điều kiện còn lại được thể hiện ở Định lý 1.5 (xem [2] trang 208).

Hình 1.1: Tập sao cấp chính xác của xấp xỉ Padé với k = 1, j = 2.

Hình 1.2: Tập sao cấp chính xác của xấp xỉ Padé với k = 2, j = 21.

Định lý 1.5 (Butcher 1964). Nếu bộ hệ số bi, ci, aij của phương pháp Runge- Kutta thỏa mãn B (p) , C (η) , D (ζ) với p ≤ η + ζ + 1 và p ≤ 2η + 2 thì phương pháp có cấp p.

Các phương pháp Gauss

Các phương pháp Gauss hay còn gọi "các phương pháp Kuntzmann-Butcher" là các phương pháp trùng khớp dựa trên cơ sở của công thức cầu phương Gauss với các hệ số c1, ..., cs là nghiệm của đa thức trực giao Legendre bậc s

ds dxs

(cid:0)xs(x − 1)s(cid:1) .

(b) p = 4

√

−

1 2

1 4

3 6

1 4

√

(a) p = 2 1 2 1

3 6 √ 3 6

1 2 1 2

3 6

1 4

1 2

1 4 1 2

Bảng 1.5: Các phương pháp Gauss cấp 2 và cấp 4.

Một số phương pháp Gauss được thể hiện trong Bảng 1.5.

Định lý 1.6 (Butcher 1964, Ehle 1968). Phương pháp Gauss s nấc có cấp chính xác p = 2s. Hàm ổn định là xấp xỉ Padé (s, s) và phương pháp ổn định A.

Các phương pháp Radau IA và Radau IIA

I :

Chứng minh. Chi tiết xem [2] và [3].

(cid:0)xs(x − 1)s−1(cid:1) , (cid:0)Radau trái(cid:1) (1.19)

II :

(cid:0)xs−1(x − 1)s(cid:1) , (cid:0)Radau phải(cid:1) (1.20)

III :

(Lobatto)

(cid:0)xs−1(x − 1)s−1(cid:1) . (1.21) Butcher (1964) đã giới thiệu các phương pháp Runge-Kutta dựa trên cơ sở của các công thức cầu phương Radau và Lobatto. Ông ấy gọi chúng theo ba loại I, II hoặc III tùy thuộc vào c1, c2, ..., cs là nghiệm của ds−1 dxs−1 ds−1 dxs−1 ds−2 dxs−2

Các trọng số bi được chọn để công thức cầu phương thỏa mãn B (s), trong đó B (2s − 1) là cơ sở Radau, B (2s − 2) là cơ sở Lobatto. Ehle (1969) đã theo đuổi ý tưởng của Butcher và xây dựng các phương pháp loại I, II và III với các tính chất ổn định vượt trội. Một cách độc lập, Axelsson (1969) đã tìm ra các phương pháp Radau IIA và chứng minh tính ổn định A của chúng.

Phương pháp Radau IA s nấc thuộc loại I khi các hệ số aij (i, j = 1, ..., s) được định nghĩa bởi điều kiện D (s). Đây là cách duy nhất bởi vì các ci là phân biệt

(a) p = 1

√

−

1 1

√

6 −

2 3

√

6 +

(b) p = 3 1 4 1 4 1 4

1 4 5 12 3 4

−1 − 18 88 + 7 360 88 + 43 360

−1 + 18 88 − 43 360 88 − 7 360

√

16 + 36

1 9 1 9 1 9 1 9

16 − 36

Bảng 1.6: Một số phương pháp Radau IA

và bi (cid:54)= 0. Bảng 1.6 trình bày các phương pháp đầu tiên có đặc điểm này. Công thức loại II của Ehle thu được bằng cách áp dụng điều kiện C (s).

(a) p = 3

(b) p = 5

√

4 −

1 3

88 − 7 360

√

4 +

269 − 169 1800 88 + 7 360

269 + 169 1800 √

√

5 12 3 4 3 4

− 1 12 1 4 1 4

√

16 − 36 16 − 36

16 + 36 16 + 36

−2 + 3 225 −2 − 3 225 1 9 1 9

Bảng 1.7: Các phương pháp Radau IIA cấp 3 và cấp 5.

Theo Định lý II.7.7 trong [2], các hệ số của phương pháp loại II được xây dựng dựa trên cơ sở là nghiệm của (1.20). Ta gọi chúng là phương pháp Radau IIA. Với s = 1 chúng ta thu được công thức Euler ẩn (xem Bảng 1.3). Các ví dụ về phương pháp Radau IIA được đưa ra trong Bảng 1.7.

Định lý 1.7. Phương pháp Radau IA s nấc và phương pháp Radau IIA s nấc đều có cấp chính xác p = 2s − 1. Hàm xấp xỉ của chúng là xấp xỉ Padé (s − 1, s). Cả hai phương pháp đều ổn định A.

Các phương pháp Lobatto IIIA, IIIB và IIIC

Chứng minh. Chi tiết xem trong [3].

i = 1, 2, ..., s

Với tất cả các công thức loại III thì ci là nghiệm của đa thức cho bởi (1.21) và các trọng số bi thỏa mãn B (2s − 2). Các hệ số aij được định nghĩa bởi C (s) cho ta các công thức Lobatto IIIA. Do đó nó là một phương pháp trùng khớp. Với các phương pháp Lobatto IIIB chúng ta áp dụng D (s). Cuối cùng, để có công thức Lobatto IIIC chúng ta đặt

ai1 = b1

(1.22)

và xác định các aij còn lại bởi C (s − 1). Ehle (1969) đã giới thiệu hai lớp phương pháp đầu tiên, và trình bày các phương pháp IIIC với s ≤ 3. Định nghĩa chung của các phương pháp IIIC được đưa ra bởi Chipman (1971) và Axelsson (1972). Các ví dụ được đưa ra ở các Bảng 1.8 - 1.10.

Định lý 1.8. Các phương pháp Lobatto IIIA, IIIB và IIIC có cấp chính xác p = 2s − 2. Hàm ổn định của các phương pháp Lobatto IIIA và Lobatto IIIB là xấp xỉ Padé (s − 1, s − 1). Với các phương pháp Lobatto IIIC, hàm ổn định là xấp xỉ Padé (s − 2, s). Và tất cả chúng đều ổn định A.

0 1 2

(a) p = 2 0 0 0 1 1 2 2 1 1 2 2

(b) p = 4 0 5 24 1 6 1 6

0 1 3 2 3 2 3

0 −1 24 1 6 1 6

Bảng 1.8: Các phương pháp Lobatto IIIA cấp 2 và cấp 4.

1 2

(a) p = 2 1 2 1 2 1 2

1 2

(b) p = 4 −1 1 6 6 1 1 3 6 5 1 6 6 1 2 3 6

1 6

Bảng 1.9: Các phương pháp Lobatto IIIB cấp 2 và cấp 4.

(a) p = 2

1 2

1 2 1 2 1 2

− 1 2 1 2 1 2

(b) p = 4 − 1 3 5 12 2 3 2 3

1 6 1 6 1 6 1 6

1 6 − 1 12 1 6 1 6

Bảng 1.10: Các phương pháp Lobatto IIIC cấp 2 và cấp 4.

Chứng minh. Chi tiết xem Phần IV.5 trong [3].

Tóm tắt về việc xây dựng các phương pháp Runge-Kutta ẩn được thể hiện

trong Bảng 1.11.

Phương pháp

Bộ điều kiện đơn giản hóa

Cấp chính xác Hàm ổn định

Gauss

B (2s)

C (s)

D (s)

(s, s) Padé

Radau IA

B (2s − 1) C (s − 1) D (s)

2s − 1

(s − 1, s) Padé

Radau IIA

B (2s − 1) C (s)

D (s − 1)

2s − 1

(s − 1, s) Padé

Lobatto IIIA B (2s − 2) C (s)

D (s − 2)

2s − 2

(s − 1, s − 1) Padé

Lobatto IIIB B (2s − 2) C (s − 2) D (s)

2s − 2

(s − 1, s − 1) Padé

Lobatto IIIC B (2s − 2) C (s − 1) D (s − 1)

2s − 2

(s − 2, s) Padé

Bảng 1.11: Thể hiện đầy đủ của các phương pháp Runge-Kutta ẩn.

W -biến đổi

√

dk (cid:16)

xk(x − 1)k(cid:17)

Ta kí hiệu

√

k (cid:88)

(−1)j+k

2k + 1

Pk (x) =

(cid:19)

2k + 1 k!

dxk

j=0

(cid:18) k j (cid:19) (cid:18) j + k j

(1.23)

1 (cid:90)

P 2

là đa thức Legendre đã được chuẩn hóa

k (x) dx = 1.

(1.24)

x (cid:90)

Các đa thức thỏa mãn

P0 (t) dt = ξ1P1 (x) +

P0 (x) ;

1 2

0 x (cid:90)

k = 1, 2...

(1.25a)

Pk (t) dt = ξk+1Pk+1 (x) − ξkPk−1 (x)

(1.25b)

với

ξk =

√ 2

4k2 − 1 Định lý 1.9 (xem [3]). Cho W định nghĩa bởi

i = 1, ..., s

j = 1, ..., s,

(1.26)

wij = Pj−1 (ci)

(1.27)

và A là ma trận hệ số của phương pháp Gauss có cấp p = 2s thì

1/2 −ξ1 ξ1

W−1AW =

 

:= XG.

0 −ξ2 . . . ξ2 . . .

(1.28)

ξs−1

          . . . 0 −ξs−1

i = 1, ..., s j = 1, ..., η + 1.

wij = Pj−1 (ci)

Bổ đề 1.2 (xem [3]). Cho A là ma trận hệ số của một phương pháp Runge-Kutta ẩn và W là ma trận không suy biến với

Khi đó điều kiện C (η) (với η ≤ s − 1) tương đương với điều kiện η cột đầu tiên của W−1AW giống với η cột đầu tiên của XG trong (1.28).

i = 1, ..., s j = 1, ..., ζ + 1,

wij = Pj−1 (ci)

Bổ đề 1.3 (xem [3]). Cho W là ma trận không suy biến với

và B = diag (b1, ..., bs) với bi (cid:54)= 0. Khi đó điều kiện D (ζ) (với ζ ≤ s − 1) tương đương với điều kiện ζ hàng đầu tiên của ma trận (cid:0)WT B(cid:1) A(cid:0)WT B(cid:1)−1 giống với ζ hàng đầu tiên của XG trong (1.28). (Nếu B suy biến thì ta vẫn có (1.29) như bên dưới).

1/2 −ξ1 0 ξ1 . . .

WT B.

WT BA =

 

(1.29)

∗ ∗

           

0 . . . . . .

−ξζ . . . . . .

∗ ∗

. . . . . . −ξζ−1 ξζ−1 . . . . . .

Hai ma trận biến đổi ở trong các Bổ đề 1.2 và Bổ đề 1.3 ở trên có thể bằng

WT B = W−1 hoặc WT BW = I.

nhau, tức là

(1.30)

Bổ đề 1.4 (xem [3]). Với bất kì công thức cầu phương mà có cấp p ≥ 2s − 2 thì ma trận

W = (Pj−1 (ci))i,j=1,...,s

(1.31)

i=1

thỏa mãn (1.30).

i = 1, ..., s,

j = 1, ..., max (η, ζ) + 1

T (η, ζ)

Định nghĩa 1.7. Với η, ζ là các số nguyên từ 0 đến s − 1. Ta nói rằng ma trận W cấp s × s thỏa mãn T (η, ζ) với công thức cầu phương có bộ hệ số (bi, ci)s nếu

a) W là ma trận không suy biến b) wij = Pj−1 (ci) (cid:18) I

c) WT BW =

0 0 R

(cid:19)

  

I là ma trận đơn vị cấp (ζ + 1) × (ζ + 1), R là ma trận tùy ý cấp (s − ζ − 1) ×

(s − ζ − 1) .

i=1 thì một phương pháp Runge-Kutta dựa trên (bi, ci)s

Định lý 1.10 (xem [3]). Với W thỏa mãn T (η, ζ) cho công thức cầu phương (bi, ci)s i=1, với ma trận X = W−1AW ta có

a) η cột đầu tiên của X là của XG ⇔ C (η);

b) ζ hàng đầu tiên của X là của XG ⇔ D (ζ).

Bổ đề 1.5 (xem [3]). Nếu công thức cầu phương có các nút ci, tất cả các trọng số bi > 0 và nếu nó có cấp chính xác p thỏa mãn p ≥ 2η + 1, p ≥ 2ζ + 1 thì ma trận

W = (pj−1 (ci))i,j=1,...,s

(1.32)

s (cid:88)

(cid:104)p, r(cid:105) =

bip (ci) r (ci).

i=1

có tính chất T (η, ζ) và thỏa mãn (1.30). Ở đây, pj (x) là đa thức trực giao bậc j với tích vô hướng

Nhận xét 1.2 (xem [3]). Chúng ta có thể biểu diễn hàm ổn định của phương pháp Runge-Kutta ẩn dưới dạng ma trận Runge-Kutta chuyển vị X = W−1AW. Từ (b) và (c) của tính chất T (η, ζ), ta suy ra

We1 = 111 , WT B111 = e1 ,

e1 = (1, 0, ..., 0)T .

(1.33)

R (z) = 1 + zeT

Do đó các công thức (1.10) và (1.11)trở thành

1 (I − zX)−1e1 , (cid:1)

R (z) =

(1.34)

det (cid:0)I − zX + ze1eT det (I − zX)

1.3 Áp dụng các phương pháp Runge-Kutta giải bài toán

cương

(1.35)

Để minh họa cho các phương pháp Runge-Kutta được xây dựng trong phần trên, ta áp dụng các phương pháp Runge-Kutta giải một số bài toán cương. Trước hết ta tìm hiểu xem bài toán cương là gì? Có nhiều cách định nghĩa khác nhau về bài toán cương. Tuy nhiên, quan điểm có tính thực tiễn nhất

cũng là quan điểm đầu tiên trong lịch sử về bài toán cương là do Curtiss và Hirschfelder đưa ra. Theo đó, bài toán cương là những bài toán mà các phương pháp Runge-Kutta hiển không giải quyết được. Đối với các bài toán cương, các phương pháp Runge-Kutta ẩn giải quyết tốt hơn rất nhiều so với các phương pháp Runge-Kutta hiển.

Định nghĩa 1.8. Bài toán giá trị ban đầu được gọi là bài toán cương nếu điều kiện ổn định tuyệt đối xác định bước đi h nhỏ hơn nhiều so với yêu cầu về độ chính xác.

Ngoài ra, ta có thể nhận dạng bài toán cương trong một vài trường hợp đặc

biệt. Với bài toán (1.6) ta có định nghĩa sau

|Reλi|

Định nghĩa 1.9. Đối với bài toán y(cid:48) = Ay, A ∈ Cm×m nếu Re (λi) ≤ 0 (với λi là

(cid:29) 1 thì bài toán gọi là bài

|Reλi|

max i min i

các giá trị riêng của A) với i = 1, ..., m và

toán cương.

Bài toán 1.1. Bài toán Van der Pol

0 ≤ x ≤ 3.

y(cid:48) 1 = y2 (cid:1) y2 − y1 2 = (cid:0)(cid:0)1 − y2 y(cid:48) y1 (0) = 2, y2 (0) = 0

  (1.36) (cid:1) /ε 

Ta xét bài toán với ε = 10−3. Ta lập trình chương trình thử nghiệm số để giải bài toán Van der Pol bằng các phương pháp Gauss cấp 2 và cấp 4 (xem Bảng 1.5) trong môi trường Matlab. Sau đó, ta so sánh kết quả lời giải số của chương trình trên với lời giải số của chương trình ode23s. Chương trình ode23s được xây dựng dựa trên cặp công thức Rosenbrock cấp 2 và cấp 3 (xem [6]). Kết quả các lời giải số của y1 thể hiện trong Hình 1.3. Sau 2000 bước với bước đi đều h = 0.0015, phương pháp Gauss cấp 4 cho ta lời giải khá chính xác, còn lời giải của phương pháp Gauss cấp 2 có nhiều sai số. Trong khi đó, nếu ta sử dụng chương trình ode23s thì chỉ cần 450 bước với hmin ≈ 1.4629 × 10−6, hmax ≈ 0.0431 ta đã có lời giải khá chính xác.

Bài toán 1.2. Bài toán cương Robertson

3.107y2 2

  (1.37)

y(cid:48) 1 = −0.04y1 + 104y2y3 y(cid:48) 2 = 0.04y1 − 104y2y3 − 3.107y2 2 y(cid:48) 3 = y1 (0) = 1; y2 (0) = 0; y3 (0) = 0, 0 ≤ x ≤ 40.



Hình 1.3: Lời giải của y1 khi ε = 10−3 bằng các phương pháp Gauss và chương trình ode23s.

Hình 1.4: Lời giải của y2 bằng các phương pháp Radau IIA cấp 3, cấp 5 và chương trình ode23t.

Ta lập trình chương trình trong Matlab sử dụng các phương pháp Radau IIA cấp 3 và cấp 5 (xem Bảng 1.6) với bước đi đều h = 0.01 để giải. Tiếp theo, ta so sánh với lời giải của chương trình ode23t. Chương trình ode23t là chương trình cải biên dựa trên công thức hình thang sử dụng các bước đi h thay đổi (xem [6]). Các lời giải số của y2 thể hiện trong Hình 1.4. Sau 4000 bước với

bước đi đều h = 0.01, lời giải của các phương pháp Radau IIA cấp 3 và cấp 5 có độ chính xác cao. Trong khi đó chương trình ode23t chỉ cần 44 bước với hmax ≈ 4, hmin ≈ 2.0026 × 10−4 đã cho ta lời giải có độ chính xác cao.

Bài toán 1.3. Bài toán cương OREGO

(y3 − (1 + y1) y2)

(cid:1)(cid:1) (cid:0)1 − 8.375.10−6y1 − y2

(1.38)

1 = 77.27 (cid:0)y2 + y1 y(cid:48) 1 y(cid:48) 2 = 77.27 y(cid:48) 3 = 0.161 (y1 − y3) y1 (0) = 1, y2 (0) = 2, y3 (0) = 3, 0 ≤ x ≤ 400.

  

Hình 1.5: Lời giải của y2 khi sử dụng chương trình ode23s.

1.4 Các loại chuẩn

Chuẩn vectơ

Nhận thấy rằng đây là một bài toán rất cương, ta sử dụng chương trình ode23s để giải. Lời giải số của y2 khi sử dụng chương trình ode23s thể hiện trong Hình 1.5. Mặc dù đây là một bài toán rất cương, chương trình ode23s chỉ cần thực hiện 694 bước với hmax ≈ 9.6650, hmin ≈ 4.0532 × 10−4.

(cid:107)y(cid:107) =

y2 1 + ... + y2 n,

Định nghĩa 1.10. Xét vectơ y = (y1, ...., yn)T , chuẩn Euclid của y (modulus của y) là (cid:113) (1.39)

n (cid:88)

(cid:107)y + z(cid:107) ≤ (cid:107)y(cid:107) + (cid:107)z(cid:107) ,

≤

nó thỏa mãn tất cả các tính chất thông thường của một chuẩn, ví dụ như bất đẳng thức tam giác (1.40)

(cid:107)yi(cid:107).

i=1

(1.40) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)

Chuẩn Euclid cho bởi (1.39) không phải là duy nhất. Ngoài ra chúng ta còn có các chuẩn cho bởi các công thức (1.41) và (1.42). Hơn nữa, người ta đã chứng minh được rằng trong không gian hữu hạn chiều vai trò của các chuẩn là tương đương nhau.

(cid:107)y(cid:107) = max (|y1| , ... , |yn|) ,

(1.41)

(cid:107)y(cid:107) = |y1| + ... + |yn| .

Chuẩn ma trận

(1.42)

Cũng giống như chuẩn của vectơ, ma trận cũng có một số chuẩn và các chuẩn

là tương đương.

(cid:107)Qu(cid:107) .

Định nghĩa 1.11 (xem [2]). Cho Q là một ma trận (n cột, m hàng) và (cid:107).(cid:107) là một trong các dạng chuẩn (1.39), (1.41), (1.42). Chuẩn tương thích của ma trận Q được định nghĩa bởi

(cid:107)Qv(cid:107) (cid:107)v(cid:107)

(cid:107)Q(cid:107) = sup v(cid:54)=0

= sup (cid:107)u(cid:107)=1

(1.43)

(cid:107)Qv(cid:107) ≤ (cid:107)Q(cid:107) (cid:107)v(cid:107) với mọi v.

Theo định nghĩa, (cid:107)Q(cid:107) là số nhỏ nhất sao cho

(1.44)

Định lý tiếp theo sẽ đưa ra các cách tính (1.43).

Định lý 1.11 (xem [2]). Chuẩn của ma trận Q = (qki) (k = 1, m, i = 1, n) được tính bởi các công thức sau đây: Với chuẩn Euclid (1.39),

(cid:107)Q(cid:107) =

(cid:113) Giá trị riêng lớn nhất của QT Q ; (1.45)

Với chuẩn max (1.41),

(cid:107)Q(cid:107) = max

;

|qki|

k=1,...,m

i=1

(cid:33) (cid:32) n (cid:88) (1.46)

|qki|

(cid:107)Q(cid:107) = max i=1,...,n

k=1

Với chuẩn (1.42), (cid:33) (cid:32) m (cid:88) (1.47)

Chuẩn logarit

y(cid:48) = f (x, y) ,

Xét bài toán Cauchy

y (x0) = y0.

(1.48)

(x, y)

≤ L với (x, y) ∈ U

Định lý 1.12 (xem [2]). Nếu f (x, y) là khả vi theo y trong một tập mở lồi U và nếu

∂f ∂y

(1.49) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)

(cid:107)f (x, z) − f (x, y)(cid:107) ≤ L (cid:107)z − y(cid:107)

thì

với (x, y) và (x, z) ∈ U. (1.50)

(Chuẩn ở (1.49) phụ thuộc vào chuẩn sử dụng trong (1.50).)

Điều kiện (1.50) được gọi là điều kiện Lipschitz. Bài toán còn có thể viết ở

dạng hệ phương trình vi phân sau

y(cid:48) 1 = f1 (x, y1, ..., yn) , y1 (x0) = y10 ... y(cid:48) n = fn (x, y1, ..., yn) , yn (x0) = yn0

  (1.51) 

Giả sử tồn tại n lời giải y1 (x) , ..., yn (x) và sử dụng phương pháp của Euler (xem [2] Phần I.9) ta được

yk,i+1 = yki + (xi+1 − xi) fk (xi, y1i, ..., yni) k = 1, ..., n, i = 0, 1, 2, ...

yki được xấp xỉ bởi yk (xi) và x0 < x1 < x2.... Giả sử ta chọn được v (x) để

v(cid:48)

(1.52)

i (x) = f (xi, yi) với xi < x < xi+1.

(1.53)

m (x) = (cid:107)v (x) − y (x)(cid:107)

Với chuẩn bất kì ta xét sai số

(1.54)

như là một hàm số theo biến x và cố gắng ước lượng sự phát triển của sai số này. Vì m (x) không phải lúc nào cũng khả vi nên ta xét đạo hàm Dini được định nghĩa bởi

m (x + h) − m (x) h

D+m (x) = lim sup h→0, h>0

(1.55)

D+m (x) = lim inf h→0, h>0

m (x + h) − m (x) h

(1.56)

Định lý 1.13 (xem [2]). Nếu hàm g (x, y) là hàm liên tục và thỏa mãn điều kiện Lipschitz thì

⇒ m (x) ≤ u (x) với x0 ≤ x

D+m (x) ≤ g (x, m (x)) D+u (x) ≥ g (x, u (x)) m (x0) ≤ u (x0)

  (1.57) 

là đúng với các hàm liên tục m (x) và u (x).

µ (Q) = lim

Định nghĩa 1.12. Cho Q là một ma trận vuông, khi đó ta gọi

h→0,h>0

(cid:107)I + hQ(cid:107) − 1 h

(1.58)

là chuẩn logarit của Q.

Định lý 1.14 (xem [2]). Chuẩn logarit (1.58) thu được bằng các công thức sau:

Với chuẩn Euclid (1.39),

µ (Q) = λmax = giá trị riêng lớn nhất của

1 2

(cid:0)QT + Q(cid:1) ; (1.59)

Với chuẩn max (1.41),

µ (Q) = max

|qki|

k=1,2,...,n

i(cid:54)=k

  (cid:88) (1.60) qkk +  ;

Với chuẩn (1.42),

µ (Q) = max

|qki|

i=1,2,...,n

k(cid:54)=i

  (cid:88) (1.61) qii +  .

Q =

Ví dụ 1.4. Cho ma trận (cid:33)

(cid:32) 1 0 0 2 1 2 1 0 1

1 2

Với chuẩn Euclid, µ (Q) = giá trị riêng lớn nhất của (cid:0)QT + Q(cid:1) ≈ 2.6861. Với

chuẩn max, µ (Q) = 5. Với chuẩn (1.42), µ (Q) = 4.

Nhận xét 1.3. Với Q ∈ Cn×n, các công thức trên vẫn đúng nếu ta thay thế QT bởi Q∗ và qkk, qii bởi Reqkk, Reqii.

Định lý 1.15 (xem [2]). Giả sử ta có ước lượng

(x, η)

≤ l (x) với η ∈ [y (x) , v (x)] và

(cid:19) (1.62) (cid:18) ∂f ∂y

(cid:13)v(cid:48) (x + 0) − f (x, v (x))(cid:13) (cid:13) (cid:13) ≤ δ (x) , (cid:107)v (x0) − y (x0)(cid:107) ≤ (cid:37).

Khi đó với x > x0 ta có

x (cid:90)

(cid:107)y (x) − v (x)(cid:107) ≤ eL(x)

e−Lsδ (s) ds

 

x (cid:82)

l (s) ds.

(1.63) (cid:37) +  ,

với L (x) =

(cid:104)y − z, f (x, y) − f (x, z)(cid:105) ≤ l (x) (cid:107)y − z(cid:107)2.

Kết luận

Mệnh đề 1.3 (xem [2] trang 76). Với chuẩn Euclid, điều kiện (1.62) tương đương với điều kiện

Chương này đã trình bày một số khái niệm cơ bản và kết quả với mục đích sử dụng các kết quả đó trong các Chương 2 và Chương 3. Ta sử dụng kiến thức về sao cấp chính xác trong chứng minh Định lý 2.2, sử dụng Định nghĩa 1.12 và Định lý 1.15 về chuẩn logarit trong Phần 2.1. Ngoài ra, với các công thức Runge-Kutta ẩn (Gauss, Radau, Lobatto) đã xây dựng trong phần này, ta sẽ nghiên cứu về tính ổn định B và ổn định đại số của chúng trong Chương 3. Bên cạnh đó, ta đã thử nghiệm số giải một số bài toán cương sử dụng các phương pháp Runge-Kutta.

Chương 2

Tính co cho bài toán tuyến tính

Nội dung của chương này được tham khảo từ tài liệu [3]. Xuất phát từ điều kiện ổn định tuyệt đối (1.5) của bài toán (1.4), ta mở rộng đến khái niệm "tính co" cho bài toán tuyến tính (1.6). Tính co có nghĩa là hai lời giải số của cùng một phương pháp ngày càng gần nhau hơn. Nội dung chương này trình bày các khái niệm và định lý liên quan đến tính co khi xét bài toán tuyến tính. Đầu tiên, chúng ta nghiên cứu điều kiện co khi xét chuẩn Euclid (Định lý 2.1 của von Neumann). Sau đó, ta nghiên cứu tính co trong chuẩn (cid:107).(cid:107)∞ và (cid:107).(cid:107)1. Ngoài ra, ta còn nghiên cứu các hàm tăng trưởng sai số với bài toán tuyến tính. Hơn nữa, từ các ước lượng trong bài toán tuyến tính dừng, ta mở rộng đến các bài toán với nhiễu phi tuyến nhỏ.

Trước tiên, ta xét hàm ϕ (x) là một lời giải trơn của bài toán y(cid:48) = f (x, y).

y(cid:48) (x) = f (x, ϕ (x)) +

(x, ϕ (x)) (y (x) − ϕ (x)) + ...

Chúng ta tuyến tính hóa hàm f trong lân cận của (x, ϕ(x)) ta được

∂f ∂y

(2.1)

y(cid:48) (x) =

(x, ϕ (x)) y (x) + ... = J (x) y (x) + ... .

và đặt y (x) − ϕ (x) = y (x) ta có

∂f ∂y

(2.2)

y(cid:48) = Jy, y(cid:48) = Jy.

Với xấp xỉ ban đầu, ta xét ma trận Jacobian J (x) như một hằng số và bỏ qua sai số ta được

(2.3)

Nếu ta áp dụng phương pháp Euler hiển cho (2.3) ta có

ym+1 = R (hJ) ym

(2.4)

R (z) = 1 + z

với

(2.5)

n (cid:88)

Công thức (2.4) được nghiên cứu bằng cách biến đổi ma trận J thành dạng chuẩn tắc Jordan. Chúng ta giả sử rằng J là ma trận dạng đường chéo và có các vectơ riêng là v1, ... , vn và y0 được biểu diễn theo cơ sở này như sau

y0 =

αivi.

i=1

(2.6)

n (cid:88)

Thay vào (2.4) ta có

ym =

(R (hλi))mαivi ,

i=1

(2.7)

trong đó λi là các giá trị riêng. Rõ ràng ym bị chặn khi m → ∞, nếu với tất các giá trị riêng, số phức z = hλi nằm trong tập S = {z ∈ C : |R (z) ≤ 1|} = {z ∈ C : |z − (−1)| ≤ 1} là đường tròn có bán kính là 1 và tâm là −1.

Việc nghiên cứu sự ổn định có thể dựa trên cơ sở chéo hóa ma trận Jacobi J ≈ ∂f/∂y như trên. Nhưng với những bài toán có số chiều lớn, ma trận thực hiện sự biến đổi có thể là ma trận điều kiện xấu và điều đó sẽ làm mất tất cả các ước lượng đẹp.

Ví dụ 2.1. Rời rạc hóa phương trình hyperbolic

∂u ∂t

∂u ∂x

(2.8)

bằng phương pháp các đường thẳng dẫn đến

−1

y(cid:48) = Ay , A = λ

, λ =

> 0.

 

1 ∆x

1 −1

(2.9)       . . . . . .

Ma trận A có tất cả các giá trị riêng là −λ và dựa trên sự phân tích ổn định phổ có thể chỉ ra sự hệ 2.9 ổn định tiệm cận hội tụ nhanh đến 0. Nếu số chiều lớn, lời giải của (2.9) cũng như (2.8) không có tính chất này. Do vậy, ta quan tâm chính đến giới hạn nghiêm ngặt với các lời giải số (xem (2.4))

ym+1 = R (hA) ym

(2.10)

trong các chuẩn khác nhau của Rn (hoặc Cn). Ở đây R (z) là hàm ổn định của phương pháp được sử dụng. Từ (2.10) và công thức (1.44) suy ra

(cid:107)ym+1(cid:107) ≤ (cid:107)R (hA)(cid:107) (cid:107)ym(cid:107)

(2.11)

(cid:107)R (hA)(cid:107) ≤ 1.

và lời giải số có tính co khi

Hình 2.1: Lời giải số y98 bằng phương pháp Euler ẩn với các giá trị ban đầu khác nhau.

Ta xét bài toán (2.8) với 0 ≤ x ≤ 1, 0 ≤ t ≤ 1. Khi ∆x = 0.01 ứng với việc chia đoạn [0, 1] thành N = 100 đoạn con, ma trận A là ma trận vuông cấp N −1 = 99. Sau đây, ta sử dụng phương pháp Euler ẩn để giải (2.9) với bước h = 0.01 sử dụng các giá trị ban đầu khác nhau. Ta lập trình chương trình thử nghiệm số của phương pháp Euler ẩn trong môi trường Matlab. Kết quả của chương trình được thể hiện trong Hình 2.1. Trên hình 2.1 ta thấy hai lời giải số bằng phương

Hình 2.2: Độ lệch giữa hai lời giải số trong Hình 2.1.

pháp Euler ẩn ngày càng gần nhau hơn. Độ lệch giữa hai lời giải số giảm nhanh về 0, thể hiện trên Hình 2.2.

2.1 Chuẩn Euclid (Định lý von Neumann)

Xét chuẩn Euclid sinh ra bởi tích vô hướng ký hiệu là (cid:104). , .(cid:105). Khi đó, với y(cid:48) = Ay

(cid:107)y(cid:107)2 =

(cid:104)y, y(cid:105) = 2Re (cid:10)y, y(cid:48)(cid:11) = 2Re (cid:104)y, Ay(cid:105) .

ta có

d dx

(2.12)

Re (cid:104)y, Ay(cid:105) ≤ 0 , ∀y ∈ Cn.

Do đó, điều kiện co của bài toán y(cid:48) = Ay(cid:48) là

(2.13)

Kết quả trên có được từ Mệnh đề 1.3 với l (x) = µ2 (A)

Re (cid:104)y, Ay(cid:105) ≤ µ2 (A) (cid:107)y(cid:107)2,

(2.14)

ở đây µ2 (A) là chuẩn logarit của A được tính theo công thức (1.59).

|R (z)| .

Định lý 2.1. Cho hàm hữu tỉ R (z) bị chặn với Rez ≤ 0 và giả sử ma trận A thỏa mãn điều kiện (2.13). Với chuẩn Euclid sinh ra bởi tích vô hướng thì ta có

(cid:107)R (A)(cid:107) ≤ sup Rez≤0

(2.15)

Nhận xét 2.1. Đây là một kết quả trong trường hợp hữu hạn chiều của Jordan- von Neumann (1951). Một chứng minh ngắn gọn được đưa ra bởi Hairer, Bader và Lubich (1982). Ý tưởng chứng minh dưới đây là của M. Crouzeix.

Chứng minh.

|R (λi)| .

(cid:107)R (A)(cid:107) = (cid:107)QR (D) Q∗(cid:107) = (cid:107)R (D)(cid:107) = max i=1,...,n

a) Với ma trận A chuẩn tắc thì có thể biến đổi về dạng đường chéo bởi một ma trận Q duy nhất. Do đó A = QDQ∗ trong đó D = diag {λ1, ..., λn}. Trong trường hợp này ta có

Và (2.15) được suy ra từ (2.13) vì các giá trị riêng của A thỏa mãn Reλi ≤ 0.

A (ω) =

(A + A∗) +

(A − A∗) .

ω 2

1 2

b) Với ma trận A tổng quát, ta xét hàm ma trận

(cid:104)v, A (ω) v(cid:105) = ωRe (cid:104)v, Av(cid:105) + iIm (cid:104)v, Av(cid:105)

Ta thấy bản chất

ϕ (ω) = (cid:104)u, R (A (ω)) v(cid:105)

là A (ω) thỏa mãn (2.13) với mọi ω mà Reω ≥ 0, vì vậy các giá trị riêng của A (ω) cũng thỏa mãn Reλ (ω) ≤ 0 với Reω ≥ 0. Do đó, hàm hữu tỉ

(cid:107)R (A (iy))(cid:107) (cid:107)u(cid:107) (cid:107)v(cid:107)

(cid:104)u, R (A) v(cid:105) = ϕ (1) ≤ sup y∈R

ϕ (iy) ≤ sup y∈R

(u, v cố định) không có cực trị trong Reω ≥ 0. Sử dụng A (1) = A chúng ta có đánh giá sau

|R (z)| (cid:107)u(cid:107) (cid:107)v(cid:107)

≤ sup Rez≤0

(cid:104)u, Cv(cid:105) .

(2.16)

Bất đẳng thức cuối cùng của (2.16) suy ra từ phần (a), vì A (iy) là ma trận chuẩn tắc (tức là A (iy) A(iy)∗ = A(iy)∗A (iy)). Do đó, công thức (2.15) là hệ quả của (2.16) và thực tế (cid:107)C(cid:107) = sup (cid:107)u(cid:107)≤1, (cid:107)v(cid:107)≤1

|R (z)| ≤ 1.

Hệ quả 2.1. Nếu hàm hữu tỉ R (z) ổn định A thì lời giải số yn+1 = R (hA) yn là co trong chuẩn Euclid (tức là (cid:107)yn+1(cid:107) ≤ (cid:107)yn(cid:107)) khi (2.13) được thỏa mãn.

|R (z)| ≤ 1 nên sup Rez≤0

Chứng minh. Vì ổn định A có nghĩa là max Rez≤0

Do (2.13) được thỏa mãn, theo Định lý 2.1 ta có (cid:107)R (A)(cid:107) ≤ 1. Mặt khác, ta có (cid:107)yn+1(cid:107) ≤ (cid:107)R (hA)(cid:107) (cid:107)yn(cid:107) nên (cid:107)yn+1(cid:107) ≤ (cid:107)yn(cid:107). Vì vậy tính co được đảm bảo với các lời giải số.

|R (z)| .

Hệ quả 2.2. Nếu ma trận A thỏa mãn Re (cid:104)v, Av(cid:105) ≤ ν(cid:107)v(cid:107)2 với mọi v ∈ Cn thì

(cid:107)R (A)(cid:107) ≤ sup Rez≤ν

(2.17)

2.2 Hàm tăng trưởng sai số với bài toán tuyến tính

Chứng minh. Áp dụng Định lý 2.1 với (cid:101)R (z) = R (z + ν) và (cid:101)A = A − νI.

|R (z)| .

Dựa vào các ước lượng ở trên, ta định nghĩa

ϕR (x) := sup Rez≤x

(2.18)

|R (x + iy)| .

ϕR (x) = sup y∈R

Hàm ϕR (x) được gọi là hàm tăng trưởng sai số với bài toán tuyến tính. Nó là hàm liên tục và đơn điệu tăng. Nếu hàm R (z) được xét trên nửa mặt phẳng Rez < x, theo nguyên lý cực đại ta có

Các ví dụ:

1. Phương pháp Euler ẩn

R (z) =

ϕR (x) =

1 1 − z

∞ nếu 1 ≤ x.

(cid:26) R (x) nếu − ∞ < x < 1, (2.19)

|R (∞)| nếu x ≤ ξ0,

R (x)

R (z) =

2. Hàm ổn định của công thức θ (hay công thức Rosenbrock 1 nấc)

ϕR (x) =

1 θ

1 + (1 − θ) z 1 − θz

∞

≤ x,

(2.20)

nếu    nếu ξ0 ≤ x < 1 θ

2θ (1 − θ) −∞

(cid:40) 1 − 2θ khi 0 < θ < 1, với ξ0 = khi θ ≥ 1.

3. Xấp xỉ Padé (0, 2)

R (z) =

ϕR (x) =

1 − z −

nếu − ∞ < x ≤ 0,   (2.21) nếu 0 ≤ x < 1,

R (x) 1 1 − x ∞

z2 2

 nếu 1 ≤ x.

4. Xấp xỉ Padé (1, 2)

1 +

12x2 + 12x + 9 + 10x + 7

|R (x)| √ (cid:112) 3

z 3

R (z) =

ϕR (x) =

2 (2 − x)

1 −

nếu − ∞ < x ≤ ξ0,   nếu ξ0 ≤ x < 2,

∞

2z 3

z2 6



√

10.

nếu 2 ≤ x, (2.22)

ở đây ξ0 = −6 − 3

5. Xấp xỉ Padé (2, 2)

√

1 +

1 2x +

9 + 3x2

R (z) =

ϕR (x) =

3 − x

1 −

nếu − ∞ < x ≤ 0,   nếu 0 ≤ x < 3,

∞

z 2 z 2

z2 12 z2 12

 nếu 3 ≤ x.

(2.23)

R (z) = ez − Czp+1 + O (cid:0)zp+2(cid:1)

Định lý 2.2. Giả sử R (z) là ổn định A, xấp xỉ cấp p của ez, tức là

với C (cid:54)= 0. Nếu thêm điều kiện |R (iy)| < 1 với y (cid:54)= 0 và |R (∞)| < 1 thì ta có

a) Nếu p là số lẻ

ϕR (x) = ex + O (cid:0)xp+1(cid:1) với x → 0.

(2.24)

b) Nếu p chẵn, ta chỉ có (2.24) với (−1)p/2Cx > 0, nếu không

ϕR (x) = ex + O (cid:0)xr+1(cid:1) với x → ∞

(2.25)

p 2

. với các số hữu tỉ dương r ≤

|R (z)| = ex + O (cid:0)|z|p+1(cid:1) = ex + O (cid:0)xp+1(cid:1) với x → 0.

Chứng minh. Với x → 0 thì max {|R (x + iy)| ; y ∈ R} phải dần đến gốc tọa độ. Ta tiếp tục nhận thấy rằng nó phải nằm trong tập sao cấp chính xác A = {z ∈ C; |R (z)| > |ez|}. Nếu p là số lẻ, sao cấp chính xác gồm có p + 1 cánh gần gốc (Bổ đề 1.1), với z → ∞ các thành phần của A thỏa mãn |z| ≤ D |x| , D < ∞. Do đó

Lập luận tương tự với p chẵn và (−1)p/2Cx > 0. Trong trường hợp còn lại (p chẵn và (−1)p/2Cx < 0) max {|R (x + iy)| ; y ∈ R} đạt đến gần trục ảo. Chứng minh chi tiết đã được đưa ra bởi Hairer, Bader và Lubich (1982).

Định lý 2.3. Với mỗi xấp xỉ ổn định A của ez, hàm ϕR (x) là hàm siêu mũ (hàm trên mũ), tức là thỏa mãn ϕR (0) = 1 và

ϕR (x1) ϕR (x2) ≤ ϕR (x1 + x2)

(2.26)

|R (z)| = 1.

với mọi x1, x2 cùng dấu.

|R (z)| ≤ 1 nên sup Rez≤0

Chứng minh. Ổn định A là tương đương với max Rez≤0

S (z) = R (a − z) R (z)

Theo định nghĩa hàm tăng trưởng sai số ta có ϕR (0) = 1. Do đó ta chỉ cần chứng minh (2.26). Lấy x1 và x2 cố định (cùng ≤ 0 hoặc cùng ≥ 0) và giả sử rằng ϕR (x1 + x2) < ∞. Ý tưởng là xét hàm hữu tỉ

|R (a − z) R (z)| ≤ ϕR (x1 + x2) .

trong đó a ∈ C là tham số thỏa mãn Rea ≤ x1 + x2. Do tính ổn định A và ϕR (x1 + x2) < ∞, S (z) xác định trên miền 0 ≤ Rez ≤ x1 + x2 (hoặc x1 + x2 ≤ Rez ≤ 0) và |S (z)| ≤ ϕR (x1 + x2) trên biên. Bằng cách lấy maximum, chúng ta có tất cả các z được coi như nằm trong dải được xét

Chúng ta chọn z trên đường thẳng Rez = x2 sao cho |R (z)| đạt giá trị lớn nhất, ta chọn a trên đường thẳng Rea = x1 + x2 (nghĩa là Re (a − z) = x1) để |R (a − z)| trở thành lớn nhất (xét các giới hạn). Từ đó suy ra được (2.26).

Tính chất (2.26) có một cách giải thích thực tế thú vị. Xét lời giải số yn thu

được với các cỡ bước đi. Áp dụng lặp đi lặp lại (2.11) và Hệ quả 2.2 ta được

(cid:107)ym(cid:107) ≤

(cid:107)y0(cid:107) ≤

(cid:107)y0(cid:107) ,

(cid:107)R (hkA)(cid:107)

ϕR (hkµ)

k=0

(cid:33) (cid:33) (cid:32)m−1 (cid:89) (cid:32)m−1 (cid:89) (2.27)

ϕR (hkµ) ϕR (hk+1µ) ≤ ϕR ((hk + hk+1) µ) ,

nếu bài toán y(cid:48) = Ay thỏa mãn (2.14) với µ = µ2 (A). Với µ < 0 và phương pháp ổn định A thì tất cả ϕR (hkµ) < 1. Nếu |R (∞)| < 1 thì các ϕR (hkµ) chỉ dần đến 1 khi hk → 0. Bất đẳng thức (2.26) được viết như sau

2.3 Bài toán với nhiễu phi tuyến nhỏ

có nghĩa rằng thay hai bước liên tiếp bởi một bước lớn cỡ hk + hk+1 và tăng các ràng buộc (2.27). Do đó, sau khi kết hợp một số bước liên tiếp (nếu cần thiết), chúng ta có thể giả sử hk ≥ h > 0 với mọi k. Điều này cho thấy (cid:107)ym(cid:107) ≤ (cid:37)m (cid:107)y0(cid:107) với (cid:37) = ϕR (hµ) < 1. Do đó, với lưới bất kỳ x0, x1, ..., với xm → ∞ chúng ta có sự ổn định tiệm cận, nghĩa là (cid:107)ym(cid:107) → 0 với m → ∞.

y(cid:48) = Jy + g (x, y) ,

Các ước lượng ở trên đúng trong bài toán tuyến tính dừng y(cid:48) = Jy và có thể mở rộng đến các bài toán có nhiễu phi tuyến nhỏ, còn được gọi là bài toán bán tuyến tính

(2.28)

(cid:104)y, Jy(cid:105) ≤ µ(cid:107)y(cid:107)2,

trong đó

(cid:107)g (x, y) − g (x, z)(cid:107) ≤ L (cid:107)y − z(cid:107) ,

(2.29a)

(2.29b)

với µ là chuẩn logarit của J và L nhỏ.

d dx

(cid:107)y (x) − (cid:98)y (x)(cid:107)2 = 2 (cid:10)y(cid:48) − (cid:98)y(cid:48), y − (cid:98)y(cid:11) là công thức nhận được sau khi đưa vào (2.28) với y(cid:48) và (cid:98)y(cid:48), sử dụng Bất đẳng thức Cauchy-Schwarz và các ước lượng ở (2.29) ta được

Ở đây, với sự hiện diện của yếu tố phi tuyến, sự ổn định thu được bằng cách ước lượng khoảng cách của hai lời giải gần nhau y (x) và (cid:98)y (x). Thay vì (2.12) ta có

(cid:107)y (x) − (cid:98)y (x)(cid:107)2 ≤ 2 (µ + L) (cid:107)y (x) − (cid:98)y (x)(cid:107)2.

d dx

(2.30)

Do đó, điều kiện co của bài toán là µ + L ≤ 0.

Bây giờ ta cần xây dựng tính chất chung cho các lời giải số. Về nguyên tắc, các ước lượng có thể tiến hành cho tất cả các phương pháp ở chương này. Tuy nhiên do từ phần sau, ta sẽ tìm hiểu ra rất nhiều các tính chất đẹp của phương pháp Runge-Kutta ẩn, nên ở đây ta sẽ tập chung vào phương pháp Rosenbrock.

Định nghĩa 2.1 (xem [3] Phần IV.7). Một phương pháp Rosenbrock s nấc khi áp dụng cho bài toán y(cid:48) = f (y) được xác định bởi công thức

+ hJ

i = 1, ..., s

y0 +

ki = hf

αijkj

γijkj

i−1 (cid:80) j=1

i (cid:80) j=1

(cid:32) (cid:33)

y1 = y0 +

bjkj

(2.31)

s (cid:80) j=1

  

trong đó αij, γij, bj là các hệ số xác định và J = f (cid:48) (y0) .

Mỗi bước của phương pháp Rosenbrock bao gồm một hệ phương trình tuyến tính với các ẩn số ki và với ma trận I − hγiiJ. Điều quan tâm đặc biệt là phương pháp khi γ11 = ... = γss = γ, ta chỉ cần phân tích LU một lần tại mỗi bước.

Ví dụ 2.2. Xét phương pháp Rosenbrock 1 nấc

(I − γhJ) k1 = hf (x0, y0) y1 = y0 + k1

(2.32)

R (z) =

1 + (1 − γ) z 1 − γz

với γ > 0 là tham số tùy ý. Hàm ổn định của công thức là

1 . Áp dụng (2.32) cho bài toán (2.28) 2 y1 = R (hJ) y0 + (I − γhJ)−1hg (x0, y0) .

và phương pháp là ổn định A với γ ≥

(cid:13) ≤

R (µh) +

1 +

(cid:107)y1 − (cid:98)y1(cid:107) ≤

(cid:107)y0 − (cid:98)y0(cid:107) =

(cid:107)y0 − (cid:98)y0(cid:107)

Lh 1 − γhµ

h (µ + L) 1 − γhµ

(2.33) Từ Định lý von Neumann (Hệ quả 2.2) ta có (cid:107)R (hJ)(cid:107) ≤ ϕR (hµ) và (cid:13) (cid:13)(I − γhJ)−1(cid:13) (1 − γhµ)−1. với ϕR được xác định như (2.20). Nếu ta lấy lời giải số thứ hai là (cid:98)y1, cũng được xác định bởi (2.33). Sự chênh lệch của (cid:98)y1 với y1 có thể ước lượng (cid:18) (cid:19) (cid:18) (cid:19)

1 γ

khi ξ0 < hµ < với ξ0 được đưa ra ở (2.20). Do đó, tính co xảy ra với µ + L ≤ 0

như ở trên.

Với công thức Rosenbrock tổng quát cho ở (2.31) khi áp dụng cho bài toán

ki = hg (x0 + cih, ui) + hJy0 + hJ

(aij + γij) kj

i (cid:80) j=1

y1 = y0 +

ui = y0 +

biki

aijkj ;

i−1 (cid:80) j=1

s (cid:80) i=1

(2.28) ta được

Ta có sự tương đương giữa các hằng số trong công thức sau đây.

y1 = R (hJ) y0 + h

bi (hJ) g (x0 + cih, ui)

Định lý 2.4. Lời giải số của một phương pháp Rosenbrock áp dụng cho (2.28) có thể được viết như sau

ui = Ri (hJ) y0 + h

aij (hJ) g (x0 + cjh, uj) , i = 1, ..., s.

  (2.34)

s (cid:80) i=1 i−1 (cid:80) j=1



Trong đó R (z) là hàm ổn định, Ri (z) là hàm ổn định trong và bi (z), aij (z) là

1 γ

các hàm hữu tỉ chỉ đạt cực trị tại và thỏa mãn bi (∞) = 0, aij (∞) = 0.

Nhận xét 2.2. Với những phương pháp ẩn tuyến tính, lời giải số có thể tìm được thông qua (2.34) với các hàm hữu tỉ hợp lý. Do đó các nghiên cứu tiếp theo có thể được áp dụng tốt với các phương pháp này.

ϕbi (hµ) (cid:107) (cid:98)ui − ui(cid:107)

(cid:107)(cid:98)y1 − y1(cid:107) ≤ ϕR (hµ) (cid:107)(cid:98)y0 − y0(cid:107) + hL

Bây giờ ta lấy lời giải số thứ hai (cid:98)y0, (cid:98)ui, (cid:98)y1 (được định nghĩa bởi (2.34)) khác với y1 và áp dụng bất đẳng thức tam giác. Sử dụng Định lý von Neumann (Hệ quả 2.2) và các giả thiết (2.29) thì

ϕaij (hµ) (cid:107) (cid:98)uj − uj(cid:107) .

(cid:107) (cid:98)ui − ui(cid:107) ≤ ϕRi (hµ) (cid:107)(cid:98)y0 − y0(cid:107) + hL

(2.35)

s (cid:80) i=1 i−1 (cid:80) j=1

  

Thay liên tiếp bất đẳng thức thứ hai của (2.35) vào vế phải của bất đẳng thức đầu tiên.

Định lý 2.5. Với giả thiết (2.29), sự chênh lệch giữa hai lời giải số của (2.31) có thể ước lượng bởi

(cid:107)(cid:98)y1 − y1(cid:107) ≤ (ϕR (hµ) + chL) (cid:107)(cid:98)y0 − y0(cid:107) với ϕR (x) được xác định bởi (2.18) (R (z) là hàm ổn định của (2.31) và c là hằng số phụ thuộc liên tục vào hL và hµ nhưng không phụ thuộc vào (cid:107)J(cid:107) (mức độ cương của bài toán)).

(2.36)

(cid:107)(cid:98)y1 − y1(cid:107) ≤ (1 + hC∗) (cid:107)(cid:98)y0 − y0(cid:107) với hµ ≤ Const. Ở đây, C∗ là hằng số độc lập với tính cương của bài toán (2.28).

Ước lượng này là thể hiện số của tính co khi ϕR (hµ) + hL∗ ≤ 1. Trong Định lý 2.2 ta có giả định chắc chắn rằng ϕR (x) = 1 + x + o (x), vì vậy ước lượng về tính co vẫn đạt được với µ + L∗ ≤ 0. Trong bất kỳ trường hợp nào, ta có sự ổn định A là

y(cid:48) = −100y + sin y, 0 ≤ x ≤ 0.2.

Bài toán 2.1. Bây giờ, ta sẽ đưa ra một ví dụ cụ thể minh họa bài toán (2.28)

(cid:104)y, −100y(cid:105) = −100(cid:107)y(cid:107)2,

Bài toán trên thỏa mãn các điều kiện của (2.29a) và (2.29b) với µ = −100, L = 1, vì

(cid:107)sin y − sin z(cid:107) =

2 cos

≤ (cid:107)y − z(cid:107) .

sin

 

y + z 2

y − z 2

 (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)

Hình 2.3: Chương trình ode23s với giá trị ban đầu y (0) = 1, z (0) = 3.

Tiếp theo, ta sử dụng chương trình ode23s (xem [6]) để giải bài toán. Nghiệm số thu được từ chương trình với các giá trị ban đầu lần lượt là y (0) = 1, z (0) = 3 được thể hiện trong Hình 2.3. Ta nhận thấy rằng, khoảng cách giữa hai lời giải số ngày càng giảm.

Nhận xét 2.3. Vì các hàm hữu tỉ bi và aij trong (2.34) thỏa mãn bi (∞) = 0, aij (∞) = 0, (1 − γhJ) bi (hJ) , (1 − γhJ) aij (hJ) cùng bị chặn với J thỏa mãn

(2.29) và với hµ ≤ C < γ−1. Thay vì điều kiện thứ hai của (2.29), ta có thể yêu cầu

(2.37) (cid:13)(I − γhJ)−1h (g (x, y) − g (x, z))(cid:13) (cid:13) (cid:13) ≤ l (cid:107)y − z(cid:107)

2.4 Tính co trong (cid:107).(cid:107)∞ và (cid:107).(cid:107)1

và Định lý 2.5 vẫn đúng với hL được thay thế bởi l. Giả thiết (2.29) chính là (2.37) với l = hL/ (1 − γhµ). Tuy nhiên, trong một số trường hợp đặc biệt, l có thể nhỏ hơn nhiều so với hL. Các kỹ thuật liên quan đã được Hundsdorfer (1985), Strehmel và Weiner (1987) sử dụng để chứng minh tính co và sự hội tụ của các phương pháp ẩn tuyến tính. Ostermann (1988) đã áp dụng các ý tưởng này cho bài toán với nhiễu phi tuyến nhỏ, khi hL = O (cid:0)hε−1(cid:1) với ε nhỏ (ε (cid:28) h), nhưng l có thể bị chặn một cách độc lập bởi ε−1.

Việc nghiên cứu về tính co trong các chuẩn tổng quát được đưa ra bởi Spijker (1983, 1985). Các kỹ thuật chứng minh tương tự được tìm ra bởi Bolley và Crouzeix (1978), trong khi xử lý một bài toán liên quan.

Định lý tiếp theo đưa ra điều kiện cần của tính co cho phương trình đặc biệt (2.9) với một trong hai chuẩn (cid:107).(cid:107)∞ và (cid:107).(cid:107)1. Sau đó, điều kiện tương tự cũng được đưa ra đầy đủ với tất cả các chuẩn.

Định lý 2.6. Cho A là ma trận n chiều trong (2.9) với λ ≥ 0 cố định. Với hàm hữu tỉ R (z) thỏa mãn R (0) = 1 ta có

(cid:107)R (hA)(cid:107)∞ ≤ 1

với tất cả số chiều n = 1, 2, ... (2.38)

R(j) (x) ≥ 0

chỉ khi

với x ∈ [−λh, 0] và j = 0, 1, 2, ... (2.39)

(Phát biểu tương tự vẫn đúng nếu thay (cid:107).(cid:107)∞ trong (2.38) bởi (cid:107).(cid:107)1 ).

n−1 (cid:88)

R (A) =

R(j) (−λ)

(λN )j j!

j=0

Chứng minh. Đặt h = 1 và viết A = −λI + λN trong đó N là một ma trận lũy linh. Trong một chuẩn phù hợp, (cid:107)N (cid:107) nhỏ tùy ý và từ đó ta áp dụng khai triển Taylor và N n = 0

Điều này có nghĩa (ví dụ với n = 4)

R (−λ) λR(cid:48) (−λ)

R(cid:48)(cid:48) (−λ)

R(cid:48)(cid:48)(cid:48) (−λ)

λ2 2!

R (−λ)

λR(cid:48) (−λ)

R(cid:48)(cid:48) (−λ)

R (A) =

R (−λ)

λ3 3! λ2 2! λR(cid:48) (−λ)

 

R (−λ)

             

n−1 (cid:88)

≤ 1.

Áp dụng công thức (1.46) ta thấy (cid:107)R (A)(cid:107)∞ ≤ 1 (hoặc (cid:107)R (A)(cid:107)1 ≤ 1) là tương đương với

λj j!

j=0

(2.40) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)R(j) (−λ) (cid:12)

Nếu (2.40) đúng với mọi n ≥ 1 thì chuỗi

R(j) (−λ)

λj j!

j≥0

(cid:88) (2.41)

hội tụ tuyệt đối và do đó

≤

≤ 1

1 = R (0) =

R(j) (−λ)

λj j!

j≥0

(cid:88) (cid:88) (cid:12) (cid:12) (cid:12)R(j) (−λ) (cid:12) (cid:12) (cid:12)

có nghĩa R(j) (−λ) ≥ 0 với mọi j ≥ 0. Theo khai triển Taylor

R(j) (x) =

R(k) (−λ)

(x + λ)k−j (k − j)!

k≥j

(cid:88)

bao gồm x ≥ −λ không âm, từ đó ta có (2.39).

Định lý tiếp theo sẽ chỉ ra điều kiện (2.39) là đủ cho tính co với chuẩn tùy ý. Nó dễ dàng được áp dụng cho hệ (2.9) vì ma trận của nó thỏa mãn (cid:107)A + λI(cid:107)∞ = λ.

(cid:107)A + λI(cid:107) ≤ λ

Định lý 2.7. Xét một chuẩn tùy ý và cho A với λ ≥ 0 thỏa mãn

(2.42)

R(j) (x) ≥ 0 với x ∈ [−(cid:37), 0] và j = 0, 1, 2, ... .

Nếu hàm ổn định của một phương pháp thỏa mãn R (0) = 1 và

(2.43)

thì lời giải số có tính co (tức là (cid:107)R (hA)(cid:107) ≤ 1) khi hλ ≤ (cid:37).

Chứng minh. Đặt h = 1. Do 0 ≤ λ ≤ (cid:37) ta có R(j) (−λ) ≥ 0 với mọi j. Hàm

R (z) =

R(j) (−λ)

(λ + z)j j!

j≥0

(cid:88) (2.44)

thỏa mãn |R (z)| ≤ R (−λ + r) với tất cả số phức z nằm trong hình tròn |z + λ| ≤ r. Tính chất này và (2.43) có nghĩa là không có cực trị nào của R (z) có thể nằm trong |z + λ| ≤ λ. Do đó, bán kính hội tụ của (2.44) lớn hơn hẳn λ. Do đó, từ (2.42) ta có

R (A) =

R(j) (−λ)

(A + λI)j j!

j≥0

(cid:88) (2.45)

2.5 Hệ số ngưỡng

Áp dụng bất đẳng thức tam giác cho (2.45) ta được điều cần phải chứng minh.

Số (cid:37) lớn nhất thỏa mãn (2.43) được gọi là hệ số ngưỡng của R (z).

R(j) (x) =

j = 0, 1, 2, ...,

j! (1 − x)j+1 ,

Ví dụ 2.3. Phương pháp Euler ẩn với

thỏa mãn (2.43) với mọi (cid:37) > 0. Nó có một hệ số ngưỡng là (cid:37) = ∞.

+ ... +

Rk0 (z) = 1 + z +

z2 2

zk k!

Ví dụ 2.4 (Hệ số ngưỡng với xấp xỉ Padé). Các đạo hàm của các hàm đa thức

được tính toán một cách dễ dàng. Khi xét điều kiện (2.43) cho Rk0 (z), ta chỉ cần quan tâm đến thành phần 1 + z, do đó (cid:37) = 1 với mọi k.

Với xấp xỉ Padé Rk1 (z) chỉ có một điểm kì dị đơn, do đó chúng ta có thể viết

+ đa thức của z,

Rk1 (z) =

a 1 − bz

dưới dạng

chỉ có đạo hàm hữu hạn và có thể đổi dấu (xem Ví dụ 2.3). Các giá trị số thu được thể hiện trong Bảng 2.1.

Hàm Rk2 không có điểm kì dị thực (xem Phần IV.4 trong [2]). Nhưng tính chất |R (z)| ≤ R (−(cid:37) + r) với |z + (cid:37)| ≤ r (xem chứng minh của Định lý 2.6) có nghĩa là max |R (z)| nằm trên đường tròn tâm −(cid:37) và bán kính r và nằm trên trục thực. Do đó, khi r tăng, điểm kì dị đầu tiên gặp đường tròn này phải là thực và

ở bên phải của −(cid:37). Ở đây, điều này là không thể vì xấp xỉ Rk2 (z) không bao giờ thỏa mãn tính chất (2.43). Điều này thể hiện bởi các dấu * trong Bảng 2.1.

j = 0 __

2.196

2.35

2.477

2.586

2.682

j = 1 ∞

j = 2

0.584

1.195

1.703

2.208

2.710

3.212

3.713

j = 3

j = 4

j = 5

0.353

0.770

1.081

1.424

1.794

2.185

2.590

Bảng 2.1: Hệ số ngưỡng của các xấp xỉ Padé.

Tất cả các giá trị của Bảng 2.1 được tính toán sử dụng sự phân tích của R (z) và được trích dẫn từ Kraaijevanger (1986) và van de Griend và Kraaijevanger (1986). Xem xét kĩ trong Bảng 2.1, các phương pháp với hệ số ngưỡng lớn nhất

không ổn định A. Một ngoại lệ là phương pháp Euler ẩn (k = 0, j = 1) có hệ số ngưỡng là (cid:37) = ∞.

Một nghiên cứu chi tiết về các hàm thỏa mãn (2.43) đã được đưa ra bởi S.Bernstein (1914) và được nối tiếp bởi F.Hausdorff (1921). Những hàm thỏa mãn (2.43) được gọi là hàm đơn điệu tuyệt đối trong [−(cid:37), 0]. Sau đó, S.Bernstein (1928) đã đưa ra các đặc điểm tiếp theo của các hàm đơn điệu tuyệt đối trong (−∞, 0].

∞ (cid:90)

extdα (t) ,

R (x) =

Định lý 2.8 (Bernstein 1928). Điều kiện cần và đủ để R (x) là hàm đơn điệu tuyệt đối trong (−∞, 0] là

(2.46)

trong đó α (t) là hàm bị chặn, không giảm và tích phân hội tụ với −∞ < x ≤ 0.

Đây là một kết quả chìa khóa đi đến hai định lý tiếp theo. Về cơ bản các định lý này không dễ chứng minh. Ta có thể tham khảo chứng minh của S.Bernstein (1928). Từ kết quả này ta có ngay "trường hợp giới hạn λ → ∞" của Định lý 2.7 với chuẩn tùy ý.

(cid:107)R (A)(cid:107) ≤ 1.

Định lý 2.9. Cho hàm R (x) là đơn điệu tuyệt đối trong (−∞, 0], R (0) = 1 và A là ma trận với chuẩn logarit không dương µ (A) ≤ 0 thì

∞ (cid:90)

dα (t) = R (0) = 1

≤

eAtdα (t)

(cid:107)R (A)(cid:107) =

(cid:13) ≤ 1 với x ≥ 0. Từ đó

(cid:13) dα (t) ≤ (cid:13)eAt(cid:13) (cid:13)

Chứng minh. Từ Định lý 1.15 ta có lời giải y (x) = eAxy0 của y(cid:48) = Ay thỏa mãn (cid:13)eAx(cid:13) (cid:107)y (x)(cid:107) ≤ (cid:107)y0(cid:107), do đó cũng có (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)

vì α (t) là một hàm không giảm.

Kết quả sau chứng tỏ rằng không có phương pháp Runge-Kutta nào cấp p > 1

có hàm ổn định R (x) là đơn điệu tuyệt đối trong (−∞, 0].

+ O (cid:0)x3(cid:1)

R (x) = 1 + x +

Định lý 2.10. Nếu hàm R (x) là đơn điệu tuyệt đối trong (−∞, 0] và

x2 2

với x → 0,

∞ (cid:82)

extdα (t), từ đó

thì R (x) = ex.

∞ (cid:90)

R(j) (0) =

tjdα (t).

Chứng minh. (Bolley và Crouzeix 1978) Theo (2.46) có R (x) =

∞ (cid:90)

(1 − t)2dα (t) = 0.

Vì R (0) = R(cid:48) (0) = R(cid:48)(cid:48) (0) = 1 nên

Kết luận

Do đó, hàm α (t) phải là hàm Heaviside (α (t) = 0 với t ≤ 1 và α (t) = 1 với t > 1). Thay vào (2.46) ta được R (x) = ex.

Chương này đã trình bày tính co khi xét bài toán tuyến tính. Trong suốt chương này, ta xét bài toán đặc biệt (2.9). Ban đầu, khi xét với chuẩn Euclid ta có các điều kiện liên quan đến tính co trong công thức (2.13), Định lý 2.1, Hệ quả 2.1, Hệ quả 2.2. Tiếp theo, ta cũng chỉ ra điều kiện cần của tính co cho bài toán (2.9) với hai chuẩn (cid:107).(cid:107)∞ và (cid:107).(cid:107)1 là Định lý 2.6. Với chuẩn tổng quát ta có Định lý 2.7. Ngoài ra, ta cũng đã mở rộng điều kiện co cho bài toán với nhiễu phi tuyến nhỏ. Hơn nữa, ta cũng đưa ra các định nghĩa về hàm tăng trưởng sai số với bài toán tuyến tính ϕR (x). Bên cạnh đó, ta đã thử nghiệm số giải bài toán (2.9) và một bài toán bán phi tuyến.

Chương 3

Tính ổn định B và tính co

3.1 Điều kiện Lipschitz một phía

Trong chương này, ta sẽ nghiên cứu sự ổn định của các phương pháp Runge- Kutta cho hệ phi tuyến tổng quát. Nội dung của chương này được tham khảo trong tài liệu [3]. Phần lớn các ý tưởng được lấy từ bài báo của Dahlquist (1963) với sự tổng quát khi nghiên cứu tính ổn định A của các bài toán phi tuyến. Nghiên cứu đầy đủ của Dahlquist về một lớp các hệ phi tuyến cuối cùng cũng đã thành công sau 12 năm. Trong buổi nói chuyện của Dahlquist tại hội nghị Dundee vào tháng 7 năm 1975, ông đã xét các phương trình vi phân thỏa mãn điều kiện Lipschitz một phía, và trình bày một vài kết quả đầu tiên về phương pháp đa bước. J.C Butcher (1975) sau khi tham dự hội nghị đã phát triển ý tưởng đó với các phương pháp Runge-Kutta ẩn và cho ra đời khái niệm ổn định B.

y(cid:48) = f (x, y)

Xét phương trình vi phân phi tuyến

(3.1)

(cid:104)f (x, y) − f (x, z) , y − z(cid:105) ≤ ν(cid:107)y − z(cid:107)2,

với chuẩn Euclid thì điều kiện Lipschitz một phía là

(3.2)

trong đó số ν là hằng số Lipschitz một phía của f. Giả sử rằng (3.1) có nghiệm, khi đó ta có kết quả sau đây:

(cid:107)y (x) − z (x)(cid:107) ≤ (cid:107)y (x0) − z (x0)(cid:107) eν(x−x0) với x ≥ x0.

Bổ đề 3.1. Cho hàm f (x, y) liên tục và thỏa mãn (3.2). Khi đó, với bất kỳ hai lời giải y (x) và z (x) của (3.1) ta có

m(cid:48) (x) = 2 (cid:104)f (x, y (x)) − f (x, z (x)) , y (x) − z (x)(cid:105) ≤ 2νm (x)

Chứng minh. Lấy đạo hàm m (x) = (cid:107)y (x) − z (x)(cid:107)2 ta được

m (x) ≤ m (x0) e2ν(x−x0) với x ≥ x0,

Bất đẳng thức đạo hàm này dẫn đến (xem Định lý 1.13)

từ đó suy ra m (x) ≤ (cid:107)y (x0) − z (x0)(cid:107)2e2ν(x−x0) với x ≥ x0.

≤ ν (Mệnh đề 1.3). Nếu f là hàm khả vi liên tục, thì Bổ đề

a) Trong tập lồi và mở, điều kiện (3.2) tương đương với điều (cid:19) kiện µ Nhận xét 3.1. (cid:18)∂f ∂y

3.1 trở thành một trường hợp đặc biệt của Định lý 1.15.

Re (cid:104)f (x, y) − f (x, z) , y − z(cid:105) ≤ ν(cid:107)y − z(cid:107)2,

y, z ∈ Cn,

b) Nếu y và f nhận giá trị phức thì điều kiện (3.2) được thay bởi điều kiện

(3.3)

3.2 Ổn định B và ổn định đại số

và Bổ đề 3.1 vẫn đúng.

s (cid:88)

Bất cứ khi nào ν ≤ 0 ở (3.2), khoảng cách giữa hai lời giải bất kỳ của (3.1) là hàm không tăng của x. Tính chất tương tự cũng được mong đợi từ các lời giải số. Ở đây, ta xét các phương pháp Runge-Kutta ẩn

y1 = y0 + h

bif (x0 + cih, gi)

i=1 s (cid:88)

i = 1, 2, ..., s

(3.4a)

gi = y0 + h

aijf (x0 + cjh, gj)

j=1

(3.4b)

(cid:104)f (x, y) − f (x, z) , y − z(cid:105) ≤ 0 ,

Định nghĩa 3.1 (Butcher 1975). Giả sử phương pháp Runge-Kutta thỏa mãn điều kiện co

(3.5)

(cid:107)y1 − (cid:98)y1(cid:107) ≤ (cid:107)y0 − (cid:98)y0(cid:107) với mọi h ≥ 0 đủ nhỏ.

thì lời giải số thỏa mãn

Khi đó, phương pháp được gọi là ổn định B (B-ổn định). Ở đây, y1 và (cid:98)y1 là các xấp xỉ số sau một bước xuất phát với giá trị ban đầu y0 và (cid:98)y0 tương ứng.

Rõ ràng, ổn định B kéo theo ổn định A. Điều này được thấy bằng cách áp dụng Định nghĩa 3.1 ở trên cho y(cid:48) = λy với λ ∈ C, hoặc chính xác hơn là áp dụng cho (cid:19) (cid:19) (3.6) (cid:18) α −β α β (cid:19) (cid:18) y1 y2 (cid:18) y(cid:48) 1 y(cid:48) 2

Định nghĩa 3.2. Với số nguyên dương s và c1, ..., cs là các số thực khác nhau (thường 0 < ci < 1), đa thức trùng khớp u (x) bậc s được định nghĩa bởi

i = 1, 2, ..., s

(3.7a)

u (x0) = y0 u(cid:48) (x0 + cih) = f (x0 + cih, u (x0 + cih))

(3.7b)

Lời giải số được xác định bằng công thức y1 = u (x0 + h).

Ví dụ 3.1. Với các phương pháp dựa trên cơ sở của phép cầu phương Gauss, việc chứng minh tính ổn định B là có thể (Wanner 1976). Ta ký hiệu u (x) và (cid:98)u (x) là các đa thức trùng khớp với giá trị ban đầu y0 và (cid:98)y0 và hàm khả vi m (x) = (cid:107)u (x) − (cid:98)u (x)(cid:107)2. Tại các điểm trùng khớp ξi = x0 + cih ta có m(cid:48) (ξi) = 2 (cid:104)f (ξi, u (ξi)) − f (ξi, (cid:98)u (ξi)) , u (ξi) − (cid:98)u (ξi)(cid:105) ≤ 0. Do phép cầu phương Gauss tính chính xác tích phân của đa thức m(cid:48) (x) (m(cid:48) (x)

x0+h (cid:82)

m(cid:48) (x) dx

(cid:107)y1 − (cid:98)y1(cid:107)2 = m (x0 + h) = m (x0) +

= m (x0) + h

bim(cid:48) (x0 + cih) ≤ m (x0) = (cid:107)y0 − (cid:98)y0(cid:107)2.

s (cid:80) i=1

là đa thức bậc 2s − 1) và trọng số bi > 0, ta có

Một "tiêu chuẩn đại số" cho tính ổn định B đã được chỉ ra một cách độc lập bởi Burrage và Butcher (1979) và Crouzeix (1979). Kết quả đó như sau:

Định lý 3.1. Nếu bộ hệ số của một phương pháp Runge-Kutta (3.4) thỏa mãn

i) bi ≥ 0 với i = 1, 2, ... .

i,j=1 là xác định không âm,

ii) M = (mij) = (biaij + bjaji − bibj)s

thì phương pháp là ổn định B.

Định nghĩa 3.3. Một phương pháp Runge-Kutta thỏa mãn (i) và (ii) của Định lý 3.1 gọi là ổn định đại số.

∆y0 = y0 − (cid:98)y0 , ∆y1 = y1 − (cid:98)y1 , ∆gi = gi − (cid:98)gi, ∆fi = h (f (x0 + cih, gi) − f (x0 + cih, (cid:98)gi)) ;

Chứng minh. Chứng minh Định lý 3.1. Chúng ta đưa vào các ký hiệu sau

s (cid:88)

Trừ hai công thức Runge-Kutta (3.4) tương ứng với y và (cid:98)y ta được

∆y1 = ∆y0 +

bi∆fi

i=1 s (cid:88)

(3.8a)

aij∆fj

∆gi = ∆y0 +

j=1

(3.8b)

s (cid:88)

Bình phương công thức (3.8a)

(cid:107)∆y1(cid:107)2 = (cid:107)∆y0(cid:107)2 + 2

bi (cid:104)∆fi, ∆y0(cid:105) +

bibj (cid:104)∆fi, ∆fj(cid:105)

i=1

j=1

(3.9)

aij∆fj, thay vào (3.9) ta có

s (cid:80) j=1

s (cid:88)

Từ (3.8b) ta có ∆y0 = ∆gi −

(cid:107)∆y1(cid:107)2 = (cid:107)∆y0(cid:107)2 + 2

bi (cid:104)∆fi, ∆gi(cid:105) −

mij (cid:104)∆fi, ∆fj(cid:105).

i=1

j=1

(3.10)

mij (cid:104)∆fi, ∆fj(cid:105) ≥ 0 nên ta suy ra

s (cid:80) i=1

s (cid:80) j=1

Từ (3.5) ta có (cid:104)∆fi, ∆gi(cid:105) ≤ 0, mặt khác

được phương pháp là ổn định B.

√

(cid:32)

(cid:33) 3

1 − γ

γ =

3 ± 6

1 − 2γ 1 2

γ 1 2

Bảng 3.1: Phương pháp SDIRK cấp 3.

> 0.

Ví dụ 3.2. Xét phương pháp SDIRK cho trong Bảng 3.1 dưới đây

1 2

i) b1 = b2 =

= γ −

m11 =

1 γ − 2

1 4

ii) Ta xác định ma trận M theo Định lý 3.1

(1 − 2γ) +

= −

γ −

m12 =

1 .0 − 2

1 2

(cid:18)

(1 − 2γ) −

m21 =

= γ −

m22 =

1 γ + 2 1 2 1 .0 + 2 1 γ + 2

1 2 1 γ − 2

1 4

1 4 1 . 2 1 4 1 4

(cid:19) 1 4

M =

γ −

Từ đó (cid:18) (cid:19)

(cid:19) (cid:18) 1 −1 1 1 −1 4

(cid:18) (cid:18)

λ − 2

γ −

= λ

P (λ) = |M − λI| =

γ − (cid:18)

(cid:20) (cid:18) (cid:19) 1 4 (cid:19) 1 4

γ −

−

γ −

− λ

− λ − (cid:19) (cid:18) 1 4

γ − (cid:19) 1 4

(cid:19)(cid:21) 1 4

Xét đa thức đặc trưng (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)

γ −

. Để M xác định không âm

M có các giá trị riêng là λ1 = 0, λ2 = 2

(cid:18)

(cid:19) 1 4

1 . 4

thì λi ≥ 0 (i = 1, 2) . Điều này tương đương với γ ≥

1 4

Từ (i) và (ii) thì với γ ≥ thì ta có tính chất ổn định B. Điều kiện này cũng

chính là điều kiện để phương pháp ổn định A.

Bài toán 3.1. Quay lại với bài toán (3.6) với α = −100, β = 1 ta có

, 0 ≤ t ≤ 0.1.

(cid:19) (cid:19)

(cid:18) −100 −1 −100 1 (cid:19) (cid:18) y1 y2 (cid:18) y(cid:48) 1 y(cid:48) 2

Ta áp dụng phương pháp SDIRK cấp 3 cho trong Bảng 3.1 với bước đi h = 0.01

1 ; 2

1 )T . Ta lập trình 2

và xét các giá trị ban đầu lần lượt là y (0) = (1; 1)T , z (0) = (

√

√ 3

chương trình thử nghiệm số của phương pháp SDIRK cấp 3 trong môi trường Matlab. Sự chênh lệch giữa hai lời giải số yn, zn tương ứng với các giá trị ban

3 − 6

√ 3

đầu y(0), z(0) khi γ = và γ = được thể hiện trong Hình 3.1. Hình

3 + 6 3 + 6

√

3 .

3.1 thể hiện rõ khi γ = độ chênh lệch giữa hai lời giải số giảm nhanh hơn

3 − 6

3.3 Một vài phương pháp Runge-Kutta ẩn ổn định đại số

khi γ =

• Tính dương của các trọng số của phép cầu phương (bi > 0).

• Tính chất xác định không âm của ma trận M.

Nghiên cứu chung về tính ổn định đại số được chia làm hai bước:

i=1 cấp p.

Định lý 3.2. Xét một phép cầu phương có công thức (ci, bi)s

a) Nếu p ≥ 2s − 1 thì bi > 0 với mọi i.

b) Nếu ci là nghiệm của (1.21) (phép cầu phương Lobatto) thì bi > 0 với mọi i.

√

Hình 3.1: Sự chênh lệch giữa hai lời giải số khi γ =

và γ =

3 .

3 + 6

3 − 6

Chứng minh. (Stieltjes 1884). Từ p ≥ 2s − 1, ta có tích phân của đa thức bậc 2s − 2 được tính chính xác, do đó

1 (cid:90)

dx > 0.

bi =

j(cid:54)=i

(cid:19)2 (cid:89) (3.11) (cid:18) x − cj ci − cj

Trong trường hợp của phép cầu phương Lobatto (c1 = 0, cs = 1, p = 2s − 2). Các thừa số với các chỉ số j = 1 và j = s được tính mà không cần sử dụng phép cầu phương, và lập luận tương tự ở trên ta có điều phải chứng minh.

M = BA + AT B − bbT .

Để xác minh điều kiện (ii) của Định lý 3.1, ta sử dụng WT (Phần 1.2) và xét WT M W thay vì M. Ký hiệu b = (b1, ..., bs)T , B = diag (b1, ..., bs) , A = (aij) , ta có

(3.12)

Nếu chọn W theo Bổ đề 1.5 thì WT M W = I và vì W T b = e1 = (1, 0, ..., 0)T nên điều kiện M xác định không âm tương đương với

WT M W = X + X T − e1eT 1

là xác định không âm. (3.13)

Trong đó, X = W−1AW = WT BAW như trong Định lý 1.10.

Định lý 3.3. Giả sử rằng một phương pháp Runge-Kutta với các ci phân biệt và bi > 0 thỏa mãn các điều kiện B (2s − 2) , C (s − 1) , D (s − 1) (xem (1.18)) thì phương pháp là ổn định đại số khi và chỉ khi |R (∞)| ≤ 1 (ở đây R (z) là hàm ổn định).

Chứng minh. Vì cấp của công thức cầu phương là p ≥ 2s − 2, ma trận W của Bổ đề 1.5

W = WGD , D = diag (cid:0)1, ..., 1, α−1(cid:1) ,

biP 2

(3.14)

s−1 (ci) (cid:54)= 0. Sử dụng

i,j=1 giống với (1.31) và α2 =

s (cid:80) i=1

WT BW = I ta có

D−1.

X =W−1AW =D−1W−1

GB(cid:1)−1

GBA(cid:0)WT

G AWGD = DWT

ở đây WG = (Pj−1 (ci))s

Áp dụng Bổ đề 1.2 với η = s − 1 và Bổ đề 1.3 với ζ = s − 1 ta được

X =

1/2 −ξ1 0 ξ1 . . .

 

−αξs−1 β

0 αξs−1

          . . . . . . −ξs−2 ξs−2

Ma trận này được thay vào (3.13), và ta thấy tất cả đều loại bỏ trừ β. Vì vậy điều kiện (ii) của Định lý 3.1 là tương đương với β ≥ 0. Sử dụng biểu diễn (1.35) của hàm ổn định, ta có

det

(cid:19)

|R (z)| =

det (cid:0)I − zX + ze1eT det (I − zX)

zs zs

det

− X

− X + e1eT 1 (cid:18)I z

(cid:1) (cid:18)I z (cid:19) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)

|R (∞)| =

det (cid:0)X − e1eT det X

βds−1 − α2ξ2 βds−1 + α2ξ2

s−1ds−2 s−1ds−2

Từ đó (cid:1) (3.15) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)

(trong đó dk = (k!) / (2k)!) là định thức của ma trận k chiều XG ở (1.28). Vì α2ξ2 s−1ds−2 > 0 nên |R (∞)| ≤ 1 khi và chỉ khi β ≥ 0 tương đương với phương pháp ổn định đại số.

Ta có thể so sánh các định lý với Bảng 1.11.

3.4 Ổn định AN

Định lý 3.4. Các phương pháp Gauss, Radau IA, Radau IIA, Lobatto IIIC ổn định đại số do đó cũng ổn định B.

Lý thuyết ổn định A dựa trên các phương trình tuyến tính dừng y(cid:48) = λy, trong khi ổn định B dựa trên các hệ phi tuyến tổng quát y(cid:48) = f (x, y). Câu hỏi

y(cid:48) = λ (x) y , Reλ (x) ≤ 0,

đặt ra rằng liệu có một lý thuyết ổn định nào ở giữa ổn định A và ổn định B hay không? Cách tiếp cận tự nhiên từ phương trình vô hướng, phương trình tuyến tính đến phương trình không dừng

(3.16)

trong đó λ (x) là hàm tùy ý có giá trị phức (Burrage và Butcher 1979, Scherer 1979). Kết quả trong phần này khá ngạc nhiên, tính ổn định của (3.16) với hầu hết các phương pháp Runge-Kutta là tương đương với tính ổn định B.

(1, ..., 1)T ) trở thành

Với bài toán (3.16), phương pháp Runge-Kutta (3.4) (g = (g1, ..., gs)T , 111 =

g = 111y0 + AZg , Z = diag (z1, ..., zs) ,

zj = hλ (x0 + cjh) .

(3.17)

Tính g từ (3.17) và thay vào (3.4a) ta được

y1 = K (Z) y0, K (Z) = 1 + bT Z(I − AZ)−1111

(3.18)

Định nghĩa 3.4. Một công thức Runge-Kutta được gọi là ổn định AN (AN-ổn định) nếu

|K (Z)| ≤ 1

(cid:26) với mọi Z = diag (z1, ..., zs) thỏa mãn Rezj ≤ 0

và zj = zk khi cj = ck (j, k = 1, ..., s) .

K (diag (z, z, ..., z)) = R (z)

So sánh (3.18) và (1.10) ta thấy

(3.19)

là hàm ổn định thông thường. Hơn nữa, lập luận như (3.6), ổn định B kéo theo ổn định AN. Do đó, ta có:

Định lý 3.5. Đối với các phương pháp Runge-Kutta, ta có

ổn định B ⇒ ổn định AN ⇒ ổn định A.

(f (x0, y0) + f (x1, y1)) , hàm K (Z) của

h 2

Với công thức hình thang y1 = y0 +

1 +

K (Z) =

(3.18) là

1 −

z1 2 z2 2

(3.20)

Đặt z2 = 0 và cho z1 → −∞, ta thấy phương pháp này không ổn định AN. Tổng quát hơn ta có kết quả sau đây

Định lý 3.6 (Scherer 1979). Các phương pháp Lobatto IIIA và Lobatto IIIB không ổn định AN do đó không ổn định B.

K (Z) =

Chứng minh. Như trong Mệnh đề 1.2 ta thấy

det (cid:0)I − (cid:0)A − 111bT (cid:1) Z(cid:1) det (I − AZ)

(3.21)

Theo định nghĩa, hàng đầu của ma trận A và hàng cuối của ma trận A − 111bT triệt tiêu với các phương pháp Lobatto IIIA (so sánh với chứng minh Định lý 1.8). Do đó, mẫu số của K (Z) không phụ thuộc vào z1 và tử số không phụ thuộc vào zs. Nếu chúng ta đặt z2 = ... = zs = 0 thì hàm K (Z) không bị chặn khi z1 → −∞. Điều này mâu thuẫn với ổn định AN.

Với các phương pháp Lobatto IIIB, bằng cách tương tự, cột cuối của A và cột

đầu của A − 111bT triệt tiêu.

Kết quả tiếp theo được đưa ra như đã đề cập ở trên, ổn định AN gần với ổn

định B hơn ổn định A.

Định lý 3.7 (Burrage và Butcher 1979). Giả sử

|K (Z)| ≤ 1

(cid:40)

(3.22) với mọi Z = diag (z1, ..., zs) với Rezj ≤ 0, và |zj| ≤ ε (cid:0)với ε > 0 nào đó(cid:1)

thì phương pháp là ổn định đại số (do đó cũng ổn định B).

s (cid:88)

|K (Z)|2 − 1 = 2

Chứng minh. Với ∆fi := zi∆gi và ∆y0 = 1, kết quả của (3.8) là ∆y1 = K (Z) . Chú ý rằng thực tế zi không cần là số thực, trong tính toán của chứng minh Định lý 3.1 chỉ ra rằng

biRezi|gi|2 −

mijzigizjgj ,

i=1

i,j=1 ở đây, g = (g1, ..., gs)T là lời giải của (3.17) với y0 = 1.

(3.23)

Để chứng minh bi ≥ 0, chọn zi = −ε < 0 và zj = 0 với i (cid:54)= j. Giả định (3.22) và

(3.23), tức là:

−2εbi|gi|2 − miiε2|gi|2 ≤ 0.

(3.24)

Với ε đủ bé, gi dần đến 1 và số hạng thứ hai của (3.24) là không đáng kể với bi (cid:54)= 0. Do đó bi ≥ 0.

Để xác minh điều kiện thứ hai của ổn định đại số, ta chọn zj = iεξj (ξj ∈ R).

s (cid:88)

−ε2

mijξiξj + O (cid:0)ε3(cid:1) ≤ 0.

i,j=1

Vì gi = 1 + O (ε) với ε → 0, từ (3.23) ta có

Do đó, M = (mij) là xác định không âm.

Kết hợp kết quả này với Định lý 3.1 và Định lý 3.5 ta có

Hệ quả 3.1. Với các phương pháp Runge-Kutta không suy biến (các phương pháp với tất cả các cj khác nhau) các khái niệm về ổn định AN, ổn định B và ổn định đại số là tương đương.

3.5 Các phương pháp Runge-Kutta khả quy

Một kết quả về sự tương đương giữa ổn định B và ổn định đại số với các phương pháp Runge-Kutta suy biến là rất khó chứng minh (xem Định lý 3.9 bên dưới) và đó sẽ là mục tiêu tiếp theo của chúng ta. Phần tiếp theo, ta tìm hiểu tính ổn định của các phương pháp khả quy.

Phương pháp Runge-Kutta (3.4) khi áp dụng cho tất cả các phương trình vi

phân dạng (3.1) có thể xảy ra các trường hợp sau:

i) Một vài nấc không ảnh hưởng tới lời giải số.

ii) Một số gi giống hệt nhau.

Trong cả hai trường hợp, phương pháp Runge-Kutta có thể được đơn giản hóa thành một dạng tương đương với số nấc ít hơn.

(a) DJ-khả quy 1 1 2

0 1 4 0

1 1 4 1

1 2 1 2

(b) S-khả quy 1 2 1 4 1 2

1 4 1 2

Bảng 3.2: Các phương pháp khả quy.

Minh họa cho trường hợp (i), xét phương pháp cho trong Bảng 3.2(a). Lời giải số của nó độc lập với g2 và tương đương với lời giải của công thức Euler ẩn. Với phương pháp ở Bảng 3.2(b), dễ thấy g1 = g2, khi hệ (3.4b) sở hữu một lời

giải duy nhất. Phương pháp này là tương đương với công thức trung điểm ẩn.

Trường hợp (i) ở trên còn có thể hiểu chính xác hơn như sau:

i /∈ T, j ∈ T.

Định nghĩa 3.5 (Dahlquist và Jeltsch 1979). Một phương pháp được gọi là DJ-khả quy nếu với một tập khác rỗng T ⊂ {1, 2, ..., s},

bj = 0 với j ∈ T

aij = 0 với

và (3.25)

Nếu không nó được gọi là DJ-bất khả quy.

⇒

(b) c1 A11 bT 1

(a) c1 A11 0 c2 A21 A22 bT 1

Bảng 3.3: DJ-khả quy

Điều kiện (3.25) có nghĩa là nấc j ∈ T không ảnh hưởng tới lời giải số. Điều này tốt nhất là hoán vị các nấc để các phần tử của T là các nấc cuối cùng (Cooper 1985). Sau đó bảng của phương pháp Runge-Kutta trở thành Bảng 3.3.

Một tính chất thú vị về DJ-bất khả quy và ổn định đại số của các phương

pháp Runge-Kutta được công bố bởi Dahlquist và Jeltsch (1979).

i = 1, ... , s.

bi > 0 với

Định lý 3.8. Nếu phương pháp Runge-Kutta DJ-bất khả quy và ổn định đại số thì

Chứng minh. Giả sử bj = 0 với một vài chỉ số j thì mij = 0 theo định nghĩa của M. Vì M là xác định không âm, tất cả các thành phần ở cột thứ j của ma trận M phải triệt tiêu để biaij = 0 với mọi i. Điều này có nghĩa (3.25) với T = {j|bj = 0}, mâu thuẫn với DJ-bất khả quy nên bj (cid:54)= 0 hay bj > 0.

ajk, nếu i, j ∈ Sl

aik =

k∈Sm

Định nghĩa 3.6. Một phương Runge-Kutta là S-khả quy nếu với một phân hoạch (S1, ..., Sr) nào đó của tập {1, ..., s} với r < s, ta có với mọi l và m (cid:88) (cid:88) (3.26)

Nếu không nó được gọi là S-bất khả quy. Các phương pháp vừa không phải DJ-khả quy cũng không phải S-khả quy được gọi là bất khả quy.

l ∈ Sl với l = 1, ..., r, ta xét phương pháp r nấc với các hệ số

Để hiểu điều kiện (3.26), ta giả sử rằng sau mỗi hoán vị nào đó của các nấc,

aik,

bk.

c∗ i = ci,

a∗ ij =

b∗ j =

k∈Sj

r , y∗

1, ..., g∗

1 và ta dễ dàng thấy

(cid:88) (cid:88) (3.27)

gi = g∗

y1 = y∗ 1,

l nếu i ∈ Sl,

Áp dụng phương pháp mới này vào (3.1) với g∗ rằng gi và y1 xác định bởi

(b)

⇒

1 −1

3 −2 2 −3 1 3 4 4

(a) 1 1 2 2 −1 −1 2 1 4

1 1 1 −1 −1 1 4

1 −2 −1 −2 −2 5 −3 1 1 1 4 4

Bảng 3.4: Ví dụ về phương pháp S-khả quy.

3.6 Định lý về sự tương đương giữa ổn định B và ổn định đại

số với các phương pháp S-bất khả quy

là một giải của phương pháp gốc (3.4). Với phương pháp ở Bảng 3.2(b) ta có S = {1, 2}. Một ví dụ khác về một phương pháp S-khả quy được xác định ở Bảng 3.4 (S1 = {1, 2, 3} và S2 = {4}).

Bổ đề 3.2. Cho A là ma trận hệ số của một phương pháp Runge-Kutta S-bất khả quy, khi đó tồn tại vectơ ξ ∈ Rs và η = Aξ để

(ξi − ξj) (ηi − ηj) < 0 với i (cid:54)= j.

(3.28)

ξ = 111 − εA111 với 111 = (1, ..., 1)T ,

Chứng minh. (Butcher 1982). Ý tưởng đầu tiên là đặt

η = Aξ = A111 − εA2111.

từ đó

Nếu ci (cid:54)= cj với mọi i, j thì ξi − ξj (cid:54)= 0 và ηi − ηj ngược dấu (ε đủ nhỏ), do đó ta có (3.28) đúng.

Trong việc chứng minh các trường hợp còn lại, ta xây dựng đệ quy các vectơ v0, v1, v2, ... và biểu thị bằng phân hoạch Pk của {1, 2, ..., s} xác định bởi mối quan hệ tương đương

i ∼ j ⇔ (vq)i = (vq)j , với q = 0, 1, ..., k.

(3.29)

Với một phân hoạch P nhất định của {1, 2, ..., s}, ta có không gian

X (P ) = (cid:8)v ∈ Rs; (v)i = (v)j nếu i ∼ j

đối với phân hoạch P (cid:9) .

AX (P ) (cid:54)⊂ X (P )

Với ký hiệu này, phương pháp là S-bất khả quy khi và chỉ khi

(3.30)

với tất cả các phân hoạch khác với {{1} , {2} , ..., {s}}.

vk+1 =

Ta bắt đầu với v0 = 111, P0 = {{1, 2, ..., s}} và định nghĩa

(cid:26) Avk ω nếu Avk /∈ X (Pk) nếu Avk ∈ X (Pk)

ξ = v0 − εv1 + ε2v2 − ... + (−ε)mvm

trong đó ω là một vectơ tùy ý của X (Pk) thỏa mãn Aω /∈ X (Pk). Lựa chọn để có (3.30). Sau một số hữu hạn các bước, sau m bước, ta được Pm = {{1} , {2} , ..., {s}}, vì số lượng các phần tử của Pk ngày càng tăng và tăng ngặt sau bước thứ hai. Do đó, tất cả các thành phần của vectơ

là khác nhau (với ε > 0 đủ bé) và (3.28) được thỏa mãn.

(cid:104)f (ui) − f (uj) , ui − uj(cid:105) < 0 với i (cid:54)= j.

Bổ đề 3.3. Cho u1, ..., uk và f (u1) , ..., f (uk) là các vectơ của Rn với

(cid:104)f (u) − f (v) , u − v(cid:105) ≤ 0

Khi đó, tồn tại một thác triển liên tục f : Rn → Rn thỏa mãn

với mọi u, v ∈ Rn.

< 0

γ = max i(cid:54)=j

(cid:104)f (ui) − f (uj) , ui − uj(cid:105) (cid:107)f (ui) − f (uj)(cid:107)2

Chứng minh. (Wakker 1985) Định nghĩa

f (u) =

(g (u) + u)

1 2γ

và đặt g (ui) = 2γf (ui) − ui để (cid:107)g (ui) − g (uj)(cid:107) ≤ (cid:107)ui − uj(cid:107). Một ứng dụng của Bổ đề 3.4 chỉ ra rằng tồn tại thác triển liên tục g : Rn → Rn thỏa mãn (cid:107)g (u) − g (v)(cid:107) ≤ (cid:107)u − v(cid:107) (tức là g không thể thác triển hơn nữa). Hàm

thỏa mãn được các yêu cầu.

Bổ đề 3.4 (Kirszbraun 1934). Cho u1, ..., uk và g (u1) , ..., g (uk) ∈ Rn mà

(cid:107)g (ui) − g (uj)(cid:107) ≤ (cid:107)ui − uj(cid:107) với i, j = 1, ..., k.

(3.31)

(cid:107)g (u) − g (v)(cid:107) ≤ (cid:107)u − v(cid:107) , ∀u, v ∈ Rn.

Khi đó, tồn tại một thác triển liên tục g : Rn → Rn để

(3.32)

Chứng minh. Tham khảo chứng minh của I.J.Schoenberg (1953) trong Phần IV.12 của [3].

Định lý 3.9 (Hundsdorfer và Spijker 1981). Với các phương pháp Runge-Kutta S-bất khả quy,

ổn định B ⇔ ổn định đại số.

Chứng minh. Vì Hệ quả 3.1 gần như đã bao gồm tất cả các trường hợp quan trọng trong thực tế và dễ chứng minh hơn. Định lý này có vẻ ít được quan tâm trong thực tế. Tuy nhiên, nó là một kết quả thú vị đã được nhiều người dự đoán trong nhiều năm. Để chứng minh định lý này, ta cần sử dụng đến các Bổ đề 3.2 - 3.4.

Re (cid:104)f (u) − f (v) , u − v(cid:105) ≤ 0 với mọi u, v ∈ C

Định lý 3.1 là đủ để chỉ ra rằng ổn định B và S-bất khả quy thì ổn định đại số. Với điều này, ta có s số phức z1, ..., zs thỏa mãn Rezj < 0 và |zj| < ε với ε > 0 đủ bé. Ta chỉ ra rằng tồn tại hàm liên tục f : C → C thỏa mãn

(3.33)

để các lời giải của phương pháp Runge-Kutta y1, gi và (cid:98)y1, (cid:98)gi tương ứng với y0 = 0, (cid:98)y0 = 1, h = 1, thỏa mãn

f ((cid:98)gi) − f (gi) = zi ((cid:98)gi − gi) .

(3.34)

Từ đó (cid:98)y1 − y1 = K (Z) với K (Z) xác định trong (3.18). Ổn định B thì |K (Z)| ≤ 1. Theo tính liên tục của K (Z) gần gốc tọa độ thì ta có |K (Z)| ≤ 1 với mọi zj thỏa mãn Rezj ≤ 0 và |zj| ≤ ε, để Định lý 3.7 chứng minh các kết luận.

s (cid:88)

∆gi = 1 +

aijzj∆gj

j=1

Xây dựng hàm f: ký hiệu ∆gi là lời giải của

f (gi) = tξi

(lời giải tồn tại duy nhất nếu |zj| ≤ ε và ε đủ bé). Với ξ, η như ở Bổ đề 3.2, ta định nghĩa

gi = tηi (cid:98)gi = gi + ∆gi

f ((cid:98)gi) = f (gi) + zi∆gi.

(3.35)

Re (cid:104)f (u) − f (v) , u − v(cid:105) < 0

Điều này được định nghĩa tốt khi t đủ lớn (được cố định sau đó), vì các ηi là phân biệt. Rõ ràng, gi và (cid:98)gi đại diện cho lời giải của một phương pháp Runge-Kutta với y0 = 0, (cid:98)y0 = 1 và (3.34) được thỏa mãn theo định nghĩa. Tiếp theo, ta chỉ ra rằng

nếu u (cid:54)= v (3.36)

được thỏa mãn với u, v ∈ D = {g1, ..., gs, (cid:98)g1, ..., (cid:98)gs} . Sau này, từ việc xây dựng của ξ, η, nếu u, v ∈ {g1, ..., gs}. Nếu u = gi, và v = (cid:98)gi thì đây là một hệ quả của

(cid:104)f (u) − f (v) , u − v(cid:105) = t2 (ξi − ξj) (ηi − ηj) + O (t) với t → ∞,

(3.34). Với trường hợp còn lại u = (cid:98)gi, v ∈ D\ {gi, (cid:98)gi} ta có

để (3.36) được thỏa mãn với t đủ lớn. Áp dụng Bổ đề 3.3, ta tìm được hàm liên tục f : C → C mở rộng (3.35) và thỏa mãn (3.33).

3.7 Hàm tăng trưởng sai số

Nhận xét 3.2. Butcher và Burrage (1979) phân biệt giữa ổn định BN (dựa trên hệ không dừng) và ổn định B (dựa trên hệ dừng). Vì phương trình vi phân được xây dựng trong chứng minh trên (xem (3.33)) là hệ dừng, nên hai khái niệm là tương đương với các phương pháp bất khả quy.

Tất cả các định lý ở phần trên chỉ đề cập đến tính co khi hằng số Lipschitz một phía ν ở (3.2) là 0 (xem Định nghĩa 3.1). Câu hỏi đặt ra là liệu có thể làm chặt các ước lượng khi biết rằng ν < 0 và liệu ta có thể ước lượng trong trường hợp chỉ có (3.2) với một vài ν > 0.

Định nghĩa 3.7 (Burrage và Butcher 1979). Cho ν ở (3.2) và đặt x = hν (h là cỡ bước đi). Ký hiệu ϕB (x) là số nhỏ nhất để có ước lượng

(cid:107)y1 − (cid:98)y1(cid:107) ≤ ϕB (x) (cid:107)y0 − (cid:98)y0(cid:107)

(3.37)

Re (cid:104)f (x, y) − f (x, z) , y − z(cid:105) ≤ ν(cid:107)y − z(cid:107)2.

với tất cả các bài toán thỏa mãn

(3.38)

Ta gọi ϕB (x) là hàm tăng trưởng sai số của phương pháp.

Ta xét các hàm f : R×Cn → Cn. Hàm này không mang tính tổng quát hơn (vì bất kỳ hệ nào cũng có thể viết dưới dạng thực bằng cách tách phần thực và phần ảo), nhưng nó thuận tiện hơn khi làm việc với bài toán y(cid:48) = λ (x) y ở đây λ (x) ∈ C.

y1 − (cid:98)y1 = K (Z1, ..., Zs) (y0 − (cid:98)y0)

Trong trường hợp bài toán tuyến tính không dừng y(cid:48) = A (x) y, điều kiện (3.38) trở thành µ (A (x)) ≤ ν, (trong đó µ (.) là chuẩn logarit). Đặt Zi := hA (x0 + cih), sự chênh lệch giữa hai lời giải số

ở đó

K (Z1, ..., Zs) = I + (cid:0)bT ⊗ I(cid:1) Z(I ⊗ I − (A ⊗ I) Z)−1 (111 ⊗ I)

(3.39)

và Z là ma trận đường chéo khối, với Z1, ..., Zs là các khối trên đường chéo.

Định lý 3.10. Hàm tăng trưởng sai số của một phương pháp Runge-Kutta ẩn thỏa mãn

(cid:107)K (Z1, ..., Zs)(cid:107) .

ϕB (x) =

sup µ(Z1)≤x,...,µ(Zs)≤x

(3.40)

Chứng minh. Cận trên. Sự chênh lệnh ∆y1 = y1 −(cid:98)y1 của hai lời giải Runge-Kutta thỏa mãn (3.8). Giả định (3.38) có nghĩa Re (cid:104)∆fi, ∆gi(cid:105) ≤ x(cid:107)∆gi(cid:107)2. Ta sẽ chứng minh rằng tồn tại các ma trận Zi (i = 1, ..., s) với µ (Zi) ≤ x để ∆fi = Zi∆gi. Có nghĩa là ∆y1 = K (Z1, ..., Zs) ∆y0 và như một hệ quả, vế phải của phương trình (3.40) là một cận trên của ϕB (x).

i = 2, ..., n.

Zu1 :=

Zui := xui −

f , (cid:107)g(cid:107)

(cid:104)ui, f (cid:105) u1 (cid:107)g(cid:107)

Nếu ∆gi = 0 thì ∆fi = 0 và có thể lấy một ma trận tùy ý thỏa mãn µ (Zi) ≤ x. Do đó, ta xét các vectơ f, g (g (cid:54)= 0) ∈ Cn thỏa mãn Re (cid:104)f, g(cid:105) ≤ x(cid:107)g(cid:107)2. Đặt u1 := g , cuối cùng ta được cơ sở trực giao u1, ..., un của Cn. Sau đó, ta định nghĩa (cid:107)g(cid:107) ma trận Z bởi

αiui.

n (cid:80) i=1

Ta có Zg = f, dễ thấy Re (cid:104)Zv, v(cid:105) ≤ x(cid:107)v(cid:107)2 với mọi v =

Cận dưới. Đầu tiên ta xét các phương pháp Runge-Kutta không suy biến. Với mỗi Z1, ..., Zs mà µ (Zi) ≤ x, cho A (x) là một hàm liên tục thỏa mãn hA (x0 + cih) = Zi và µ (A (x)) ≤ x với mọi x (ví dụ A (x) tuyến tính nội suy). Ta có ∆y1 = K (Z1, ..., Zs) ∆y0 suy ra ϕB (x) ≥ (cid:107)K (Z1, ..., Zs)(cid:107) với tất cả Z1, ..., Zs mà µ (Zi) ≤ x.

∆gi = ∆y0 +

aijZj∆gj có một lời giải. Chính xác, như trong Định lý 3.9 ta có

s (cid:80) j=1

Với các phương pháp suy biến việc chứng minh phức tạp hơn. Không mất tính tổng quát, ta giả sử rằng phương pháp là bất khả quy, bởi vì ϕB (x) và vế phải của phương trình (3.40) đều không thay đổi khi phương pháp được thay bởi một phương pháp tương đương. Quan tâm chính bây giờ là Bổ đề 3.3 và 3.4, với số chiều tùy ý. Xét Z1, ..., Zs với µ (Zi) ≤ x, chẳng hạn hệ tuyến tính

thể xây dựng hàm liên tục f : Cn → Cn thỏa mãn (3.38) với ν = x (ta đặt h = 1) và f (gi) − f ((cid:98)gi) = Zi (gi − (cid:98)gi) . Điều này hoàn tất chứng minh của định lý.

Với các phương pháp 1 nấc (s = 1) Định lý von Neumann (Hệ quả 2.2) là đủ để xét trong trường hợp vô hướng, giá trị phức z1 ∈ C ở phương trình (3.40). Do trong trường hợp này K (Z) = R (z), ta có

ϕB (x) = ϕR (x)

với mọi công thức 1 nấc. (3.41)

Bây giờ điều này vẫn chưa rõ ràng, tuy nhiên ta có thể hạn chế supremum ở phương trình (3.40) để có thể xét trong trường hợp vô hướng hay giá trị zi ∈ C với s ≥ 2. Điều này đòi hỏi một sự tổng quát của Định lý von Neumann với các hàm nhiều biến (Hairer và Wanner 1996). Ta sẽ quay lại câu hỏi này trong phần sau.

Định lý 3.11 (Hairer và Wanner 1996). Với các phương pháp Runge-Kutta ổn định B, hàm tăng trưởng sai số là hàm siêu mũ (hàm trên mũ), tức là ϕB (0) = 1 và

ϕB (x1) ϕB (x2) ≤ ϕB (x1 + x2)

với x1, x2 cùng dấu.

S (z) = u∗

AK (A1 − zI, ..., As − zI) vAu∗

BK (B1 + zI, ... , Bs + zI) vB,

Chứng minh. Tính chất ϕB (0) = 1 có được theo Định nghĩa 3.3. Với việc chứng minh bất đẳng thức, chúng ta xét hàm hữu tỉ

sup (cid:107)u(cid:107)=1, (cid:107)v(cid:107)=1 minh Định lý 2.3.

trong đó, các ma trận Aj, Bj thỏa mãn µ (Aj) ≤ x1+x2 và µ (Bj) ≤ 0, uA, vA, uB, vB là các vectơ tùy ý của Cn. Sử dụng tính chất µ (Aj − zI) = µ (Aj) − Rez và |u∗Cv|, thu được bất đẳng thức chính xác tương tự như chứng (cid:107)C(cid:107) =

Thực tế ϕB (x) là cận trên, cùng với ϕB (−∞) = |R (∞)| cho phép ta có những

3.8 Tính toán ϕB (x)

kết luận tương tự về sự ổn định tiệm cận của lời giải số như ở Chương 2.

Để tìm giá trị lớn nhất của (cid:107)∆y1(cid:107) dưới sự hạn chế của (3.38). Chính xác hơn,

ta xét bài toán tối ưu hóa với ràng buộc bất đẳng thức

i = 1, ..., s.

(cid:107)∆y1(cid:107)2 → max Re (cid:104)∆fi, ∆gi(cid:105) ≤ x(cid:107)gi(cid:107)2,

(3.42)

L (∆f, D) =

Ở đây, ∆f1, ..., ∆fs được coi là các giá trị độc lập trong Cn, ∆y1 và ∆gi được định nghĩa ở (3.8), và ∆y0 được xem như một tham số. Một cách tiếp cận cổ điển để giải bài toán tối ưu hóa (3.42) là xét các nhân tử Lagrange d1, ..., ds và xét

di (cid:18)(cid:18)

(∆y∗

⊗ I

s (cid:80) i=1 0, ∆f ∗)

(3.43) (cid:19) (cid:19)

α uT u W

1 (cid:107)∆y1(cid:107)2 − 2 −1 2

(cid:0)Re (cid:104)∆fi, ∆gi(cid:105) − x(cid:107)gi(cid:107)2(cid:1) (cid:19) (cid:18) ∆y0 ∆f

trong đó, ∆f = (∆f1, ..., ∆fs)T , D = diag (d1, ..., ds) và

(3.44a)

α = −1 − 2x111T D111 u = D111 − b − 2xAT D111 W = DA + AT D − bbT − 2xAT DA.

(3.44b)

(3.44c)

Định lý 3.12 (Burrage và Butcher 1980). Nếu ma trận

α + ϕ2 uT W

(cid:18) (cid:19) nửa xác định dương (3.45)

với d1 ≥ 0, ..., ds ≥ 0 thì (cid:107)∆y1(cid:107) ≤ ϕ (cid:107)∆y0(cid:107) với tất cả các bài toán thỏa mãn (3.38) với hν ≤ x. Do đó, ta có ϕB (x) ≤ ϕ.

ϕ2(cid:107)∆y0(cid:107)2 2

s (cid:88)

ta được Chứng minh. Trừ cả hai vế của (3.43) cho

1 2

i=1

(cid:0)(cid:107)∆y1(cid:107)2 − ϕ2(cid:107)∆y0(cid:107)2(cid:1) − (cid:0)Re (cid:104)∆fi, ∆gi(cid:105) − x(cid:107)∆gi(cid:107)2(cid:1) ≤ 0.

Kết quả sau đó có được do di ≥ 0 và Re (cid:104)∆fi, ∆gi(cid:105) ≤ x(cid:107)∆gi(cid:107)2.

Với sự trợ giúp của Định lý 3.12, Burrage và Butcher (1980) đã tính toán cận trên đúng của ϕB (x) cho một vài phương pháp 2 nấc. Và nhận thấy rằng, với tất cả các phương pháp 2 nấc thì ϕB (x) = ϕK (x), trong đó

|K (z1, ..., zs)| .

ϕK (x) =

sup Rez1≤x,..., Rezs≤x

(3.46)

1, ..., z0

s là các giá trị để (3.46) đạt supremum. Theo nguyên lý (cid:1) và cho ∂jK (cid:0)z0(cid:1) là

1, ... , z0 s

j = ∞). Đặt z0 = (cid:0)z0

j = x + iy0

j (y0

Có phải ϕB (x) = ϕK (x) với tất cả các phương pháp Runge-Kutta hay không? Nếu muốn kiểm tra tính hợp lệ của ϕB (x) = ϕK (x) với mỗi phương pháp Runge- Kutta, ta phải tìm các nhân tử Lagrange không âm di để (3.45) được thỏa mãn. Các bổ đề sau đây sẽ rất hữu ích cho mục đích này.

Ta ký hiệu z0 maximum ta có z0 đạo hàm của K (z1, ..., zs) theo đối số thứ j tại z0.

Bổ đề 3.5. Cho x cố định với ϕK (x) < ∞. Điều kiện (3.45) với ϕ = ϕK (x) xác định duy nhất các nhân tử Lagrange d1, ..., ds (xem phương trình (3.52) ở phần sau). Và d1, ..., ds là các số thực dương.

|K (z1, ... , zs)|2 − ϕ2 = − (1, ∆f ∗)

∆f

α + ϕ2 uT W

Chứng minh. Xét đồng nhất thức (3.43) cho trường hợp đặc biệt với ∆fj là vô hướng, ∆fj = zj∆gj do đó ∆y1 = K (z1, ... , zs). Với Rezj = x, ta có (cid:19) (cid:18) (cid:19) (cid:18) 1 (3.47)

D111 − b − 2xAT D111 + (cid:0)DA + AT D − bbT − 2xAT DA(cid:1) ∆f = 0.

Đặt ϕ := ϕK (x) và zj := z0 j (cuối cùng ta phải xét các giới hạn) vế trái của phương trình (3.47) triệt tiêu. Kết hợp với giả thiết (3.45) có nghĩa u + W∆f = 0, tức là

1, ..., z0 s

(cid:1), ta được Kết hợp hợp lý các điều đã có và sử dụng ∆f = Z0∆g , ∆g = 111 + A∆f, trong đó Z0 = diag (cid:0)z0

b.K (cid:0)z0(cid:1) .

D∆g = (cid:0)I − AT Z∗ 0

(cid:1)−1 (3.48)

Ta sẽ chứng minh tất cả các thành phần của ∆g = (I − AZ0)−1111 đều khác 0, để (3.48) xác định duy nhất các nhân tử Lagrange d1, ..., ds.

K (cid:0)z0

1 + c (cid:0)zj − z0

1, ..., zj, ..., z0 s

(cid:1)2(cid:17)(cid:17) (cid:1) + O Thác triển K (z1, ... , zs) trong chuỗi Taylor với đối số zj ta có (cid:16)(cid:0)zj − z0 (cid:1) = K (cid:0)z0(cid:1) (cid:16)

1, ..., zj, ..., z0 s

(cid:12) với Rezj ≤ (cid:12)K (cid:0)z0(cid:1)(cid:12) (cid:12) ≤ (cid:12) (cid:1)(cid:12) (cid:12)K (cid:0)z0

0 <

< ∞.

trong đó c = ∂jK (cid:0)z0(cid:1) /K (cid:0)z0(cid:1) . Vì (cid:12) Rez0 j . Ta có c > 0, do vậy ta cũng có

∂jK (cid:0)z0(cid:1) (cid:54)= 0,

∂jK (cid:0)z0(cid:1) K (z0)

(3.49)

Vi phân K (z1, ... , zs) = 1 + bT Z(I − AZ)−1111 theo zj

∂jK (cid:0)z0(cid:1) = bT (I − Z0A)−1ejeT

j (I − AZ0)−1111

(3.50)

Từ (3.49) ta có

bT (I − Z0A)−1ej (cid:54)= 0, ∆gj = eT

j (I − AZ0)−1111 (cid:54)= 0

(3.51)

để d1, ..., ds xác định duy nhất bởi (3.48). Chia thành phần thứ j của (3.48) cho ∆gj, từ (3.50)

dj = (cid:12)

K (cid:0)z0(cid:1) ∂jK (z0)

(3.52) (cid:12)bT (I − Z0A)−1ej (cid:12) (cid:12)

là số thực dương theo (3.49) và (3.51).

j là hữu hạn. Nếu z0

j = x+i∞ với một vài j nào đó, người ta phải đổi ωj = x + 1/ (zj − x), nửa mặt phẳng Rezj ≤ x thành Reωj ≤ x và ∞ thành 0.

Trong chứng minh này, ta mặc định rằng tất cả z0

Bổ đề 3.6. Nếu ma trận W của phương trình (3.44c) với d1, ..., ds như trong Bổ đề 3.5 là nửa xác định dương, thì ta có ϕB (x) = ϕK (x) .

= 0

∆f

K (x) uT α + ϕ2 W u

Chứng minh. Từ (cid:18) (cid:19) (cid:19) (cid:18) 1 (3.53)

(xem (3.47)) và vT Wv ≥ 0 với mọi v ∈ Rs, ma trận ở (3.53) nửa xác định dương. Kết quả sau đó có được từ Định lý 3.12.

• Tính ϕ = ϕK (x) của phương trình (3.46) với số hoặc với các công thức,

Từ các kết quả ở trên, với một phương pháp Runge-Kutta nhất định thì có thể kiểm tra xem khi nào ϕB (x) = ϕK (x) được thỏa mãn. Điều này có thể thực hiện bằng thuật toán sau đây:

• Tính các nhân tử Lagrange d1, ..., ds từ Bổ đề 3.5.

• Kiểm tra khi nào ma trận W của phương trình (3.44c) là nửa xác định dương. Nếu đúng W là nửa xác định dương thì ϕB (x) = ϕK (x) theo Bổ đề 3.6.

chương trình hỗ trợ.

Ví dụ 3.3. Với phương pháp Radau IIA 2 nấc p = 3 (xem Bảng 1.7)

4 5 − 2x

ϕB (x) = ϕK (x) =

nếu x ≤ ξ  

3 2

3 + 4x (cid:112)(3 − 2x) (3 + 4x − 2x2)

√

nếu ξ ≤ x < 

17(cid:1) /8.

trong đó ξ = (cid:0)9 − 3

Với phương pháp Gauss 2 nấc p = 4 (xem Bảng 1.5 )

√

1 2x +

9 + 3x2

ϕB (x) =

nếu − ∞ < x ≤ 0  

3 − x

nếu 0 ≤ x < 3. 

Với phương pháp Lobatto IIIC 2 nấc p = 2 (xem Bảng 1.10 )

1 1 − x + x2

ϕB (x) =

nếu − ∞ < x ≤ 0  

1 1 − x

nếu 0 ≤ x < 1. 

Với các phương pháp có số nấc s > 2, việc đưa ra công thức tường minh là rất khó, ta áp dụng các phương pháp số để tính toán z0 j (supremum trong phương trình (3.46)).

Kết luận

Chương này đã trình bày các khái niệm về ổn định B, ổn định đại số, ổn định AN và mối quan hệ giữa các dạng ổn định. Định lý 3.4 và Định lý 3.6 là các kết quả quan trọng đưa ra ví dụ với các phương pháp Runge-Kutta ẩn và sự ổn định đại số của chúng từ đó có kết luận về sự ổn định B. Bên cạnh đó, phần này cũng đưa ra các khái niệm về các phương pháp khả quy, các phương pháp bất khả quy. Sau đó, luận văn trình bày Định lý 3.3 về sự tương đương giữa ổn định B và ổn định đại số với các phương pháp S-bất khả quy.

Kết luận

Luận văn trình bày về tính ổn định và tính co của các phương pháp Runge- Kutta. Tính ổn định và tính co của phương pháp Runge-Kutta được trình bày lần lượt từ bài toán tuyến tính đến bài toán phi tuyến. Mối quan hệ giữa các dạng ổn định A, ổn định B, ổn định đại số, ổn định AN. Hơn nữa, chúng ta còn có các ví dụ về các phương pháp tương ứng thỏa mãn từng dạng ổn định. Đặc biệt, luận văn quan tâm đến tính ổn định của các phương pháp DJ-khả quy, S-khả quy và các phương pháp bất khả quy. Bên cạnh đó, luận văn cũng xem xét đến các hàm tăng trưởng sai số ϕR (x), ϕB (x) tương ứng khi xét bài toán tuyến tính và bài toán phi tuyến cùng với cách tính toán chúng.

Luận văn đã trình bày một số chứng minh chi tiết liên quan đến tính ổn định và tính co của các phương pháp Runge-Kutta. Các kết quả đó có thể giúp ta trong việc lựa chọn phương pháp phù hợp để giải các hệ phương trình vi phân tuyến tính và phi tuyến. Ngoài ra, luận văn có trình bày một số ví dụ, bài toán minh họa cùng với việc thử nghiệm số giải các bài toán đó.

Tài liệu tham khảo

[1] Phạm Kỳ Anh, 2008,Giải tích số, Nhà xuất bản Đại học Quốc gia Hà Nội.

[2] E. Hairer, S. P. Norsett, G. Wanner (1993), Solving Ordinary Differential

Equations I: Nonstiff Problems, Springer-Verlag, second edition.

[3] E. Hairer, S. P. Norsett, G. Wanner (1991), Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Springer-Verlag, sec- ond revised edition.

[4] J.C.Butcher (2003), Numerical Methods for Ordinary Differential Equa-

tions, The Uninversity of Auckland, New Zealand, Wiley Publishers.

[5] Linda R.Petzold, UriM.Ascher (1997), Computer Methods for Ordinary Dif-

ferential Equations and Differential-Algebraic Equations.

[6] L. Shampine, M. W. Reichelt (1997), The MATLAB ODE Suite.

[7] W. Wasow (1965),Asymptotic expansions for ordinary differential equations,

Interscience, John Wiley and Sons, New York.

Luận văn Thạc sĩ Toán học: Tính ổn định và tính co của các phương pháp Runge - kutta

Luận văn trình bày các khái niệm cơ bản về phương pháp Runge-Kutta, cách xây dựng phương pháp Runge-Kutta ẩn. Đồng thời trình bày các khái niệm và định lý liên quan đến tính co khi xét bài toán tuyến tính. Mời các bạn tham khảo!

TÍNH ỔN ĐỊNH VÀ TÍNH CO CỦA CÁC PHƯƠNG PHÁP RUNGE-KUTTA

Mục lục

Mở đầu

Bảng ký hiệu

Chương 1

Các khái niệm cơ bản

Chương 2

Tính co cho bài toán tuyến tính

Chương 3

Tính ổn định B và tính co

Kết luận

Tài liệu tham khảo

Có thể bạn quan tâm

Luận văn Thạc sĩ: Quản lý hoạt động phát triển nhận thức cho trẻ mẫu giáo ở các trường mầm non, huyện Hoàng Su Phì, tỉnh Hà Giang

Bài giảng điện tử Phương pháp tính: Số gần đúng và sai số

Bài giảng điện tử Phương pháp tính: Phương trình phi tuyến

Bài giảng điện tử Phương pháp tính: Hệ phương trình tuyến tính

Bài giảng điện tử Phương pháp tính: Đạo hàm và tích phân

Bài giảng điện tử Phương pháp tính: Nội suy và xấp xỉ hàm

Bài giảng điện tử Phương pháp tính: Phương trình vi phân thường

Bài giảng Nhóm phương pháp tính: Chương 1 - Đại học Bách khoa Hà Nội

Bài giảng Nhóm phương pháp tính: Chương 2 - Đại học Bách khoa Hà Nội

Bài giảng Nhóm phương pháp tính: Chương 3 - Đại học Bách khoa Hà Nội

Bài giảng Nhóm phương pháp tính: Chương 4 - Đại học Bách khoa Hà Nội

Bài giảng Nhóm phương pháp tính: Chương 5 - Đại học Bách khoa Hà Nội

Bài giảng Nhóm phương pháp tính: Chương 6 - Đại học Bách khoa Hà Nội

Luận văn Thạc sĩ tóm tắt: Pháp luật về điều kiện đầu tư kinh doanh trong lĩnh vực giáo dục, đào tạo ở Việt Nam

Luận văn Thạc sĩ: Sự nghiệp nghiên cứu phê bình văn học của Hoài Thanh

Luận văn Thạc sĩ: Tiểu thuyết của Đỗ Phấn từ góc nhìn sinh thái

Luận văn Thạc sĩ: Giải pháp ứng phó với nhập cư ở Liên minh Châu Âu

Luận văn Thạc sĩ: Diễn ngôn về giới nữ trong văn xuôi nữ Việt Nam đương đại (khảo sát sáng tác của Dạ Ngân, Y Ban, Lý Lan, Nguyễn Thị Thu Huệ)

Luận văn Thạc sĩ: Xây dựng mức phát thải tham chiếu rừng khu vực huyện Bảo Lâm tỉnh Lâm Đồng

Luận văn Thạc sĩ: Đặc điểm nhân vật chính trong ba tác phẩm của Franz Kafka: Lâu đài, Vụ án, Hóa thân

Tài liêu mới

Luận án Tiến sĩ: Nghiên cứu xây dựng thuật toán thích nghi và học tăng cường cấu trúc Actor - Critic điều khiển bám quỹ đạo cho robot di động đa hướng mecanum

Luận án Tiến sĩ: Cơ cấu bệnh tim mạch và chất lượng cuộc sống của người cao tuổi mắc suy tim, rung nhĩ điều trị tại Bệnh viện Thống Nhất, thành phố Hồ Chí Minh

Tóm tắt Luận án Tiến sĩ: Nghiên cứu hiện tượng nứt dăm đê sông vùng đồng bằng sông Hồng và dự báo khả năng bị nứt của một số đoạn đê

Tóm tắt Luận án Tiến sĩ: Nghiên cứu xây dựng giải pháp đảm bảo an toàn thông tin cho quá trình học liên kết dựa trên mật mã

Tóm tắt Luận án Tiến sĩ: Phát triển năng lực đánh giá công nghệ cho học sinh trong dạy học môn Công nghệ 11 ở trường trung học phổ thông

Tóm tắt Luận án Tiến sĩ: Nghiên cứu phân loại chi cầu diệp – Bulbophyllum Thouars (Orchidaceae) ở vùng Tây Nguyên bằng phương pháp hình thái và phân tử

Tóm tắt Luận án Tiến sĩ: Nghiên cứu đặc điểm phân bố và dinh dưỡng của các loài lưỡng cư ở Vườn Quốc gia Bến En và Khu bảo tồn thiên nhiên Pù Luông, tỉnh Thanh Hóa

Luận án Tiến sĩ: Tổng hợp luật dẫn và điều khiển cho một lớp tên lửa đối hải trên cơ sở ứng dụng mạng nơ ron và hệ mờ

Luận án Tiến sĩ: Nghiên cứu tổng hợp hệ điều khiển góc Pitch tua bin gió trong điều kiện có nhiễu tác động

Luận án Tiến sĩ: Nghiên cứu hóa học lipid của hai loài san hô thủy tức Millepora dichotoma và Millepora platyphylla ở Việt Nam

Luận án Tiến sĩ: Nghiên cứu kiểm soát phân phối công suất kéo trên cầu chủ động của ô tô con bằng ABS

Luận án Tiến sĩ: Ứng dụng phản ứng Domino vào tổng hợp các dẫn xuất Podophyllotoxin, Pyrimidine và đánh giá hoạt tính sinh học của các chất tổng hợp được

Luận án Tiến sĩ: Nghiên cứu thành phần hóa học và một số hoạt tính sinh học của cây chùm ngây (Moringa oleifera)

Tóm tắt Luận án Tiến sĩ: Nghiên cứu thành phần hóa học và ứng dụng ức chế ăn mòn cho thép của cao chiết xuất từ cây Lộc vừng thuộc họ Lecythidaceae

Luận án Tiến sĩ: Nghiên cứu hiện tượng nứt dăm đê sông vùng đồng bằng sông Hồng và dự báo khả năng bị nứt của một số đoạn đê

AI tóm tắt

Giới thiệu tài liệu

Đối tượng sử dụng

Từ khoá chính

Nội dung tóm tắt

Giới thiệu

Về chúng tôi

Việc làm

Quảng cáo

Liên hệ

Chính sách

Thoả thuận sử dụng

Chính sách bảo mật

Chính sách hoàn tiền

DMCA

Hỗ trợ

Hướng dẫn sử dụng

Đăng ký tài khoản VIP

093 303 0098

support@tailieu.vn

Phương thức thanh toán

Theo dõi chúng tôi

Facebook

Youtube

TikTok