Luận văn Thạc sĩ Toán học: Giải số phương trình vi phân ma trận với ràng buộc đa tạp

ĐẠI HỌC THÁI NGUYÊN

TRƢỜNG ĐẠI HỌC KHOA HỌC -------------------------------

ĐẶNG MẠNH HÙNG

GIẢI SỐ PHƢƠNG TRÌNH VI PHÂN MA TRẬN VỚI RÀNG BUỘC ĐA TẠP

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

THÁI NGUYÊN - 2019

ĐẠI HỌC THÁI NGUYÊN

TRƢỜNG ĐẠI HỌC KHOA HỌC -------------------------------

ĐẶNG MẠNH HÙNG

GIẢI SỐ PHƢƠNG TRÌNH VI PHÂN MA TRẬN VỚI RÀNG BUỘC ĐA TẠP

Chuyên ngành: Toán ứng dụng Mã số : 8 46 01 12

LUẬN VĂN THẠC SĨ TOÁN HỌC NGƯỜI HƯỚNG DẪN KHOA HỌC

TS. Nguyễn Thanh Sơn

THÁI NGUYÊN - 2019

Mục lục

Mở đầu 1

1 Kiến thức chuẩn bị 4

1.1 Một số phương pháp số giải phương trình vi phân. . . . . . . . . . . . 4

1.1.1 Phương pháp Runge-Kutta . . . . . . . . . . . . . . . . . . . . 5

1.1.2 Phương pháp Runge-Kutta phân hoạch . . . . . . . . . . . . . 6

1.1.3 Phương pháp Nystr¨om . . . . . . . . . . . . . . . . . . . . . . 7

1.2 Khái niệm đa tạp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.2.1 Đa tạp khả vi . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.2.2 Đa tạp con . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.2.3 Vectơ tiếp xúc, không gian tiếp xúc . . . . . . . . . . . . . . . 12

1.3 Đa tạp Riemann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.3.1 Khái niệm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.3.2 Khoảng cách . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Bảo toàn của tích phân và phương pháp giải 16

2.1 Khái niệm bất biến trong tích phân . . . . . . . . . . . . . . . . . . . . 16

2.2 Bất biến bậc hai . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3 Phương pháp chiếu . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.4 Tích phân phương trình vi phân trên đa tạp Stiefel . . . . . . . . . . . 26

2.4.1 Sơ lược về đa tạp Stiefel . . . . . . . . . . . . . . . . . . . . . 26

2.4.2 Cung trắc địa trên đạ tạp Stiefel . . . . . . . . . . . . . . . . . 27

i

2.4.3 Di chuyển song song . . . . . . . . . . . . . . . . . . . . . . . 29

2.4.4 Tích phân số . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.4.5 Phương pháp kiểu Runge-Kutta bậc hai . . . . . . . . . . . . . 32

2.4.6 Phương pháp giải số trên mặt cầu đơn vị . . . . . . . . . . . . 34

2.5 Phương trình vi phân trên nhóm Lie . . . . . . . . . . . . . . . . . . . 36

2.5.1 Sơ lược về nhóm Lie . . . . . . . . . . . . . . . . . . . . . . . 36

2.5.2 Phương trình vi phân trên nhóm ma trận Lie . . . . . . . . . . 38

2.5.3 Phương pháp Crouch-Grossmann . . . . . . . . . . . . . . . . 38

2.6 Ví dụ số . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Kết luận 45

Tài liệu tham khảo 46

ii

Bảng ký hiệu

⊗

C∞(M)

tích tenxơ của các không gian vectơ

Hx = ∇xH gradiant theo biến số x

tập tất cả các hàm trơn trên M

PM

sym(A)

hình chiếu trên M

skew(A)

phần đối xứng của ma trận vuông A

phần phản đối xứng của ma trận vuông A

ΠTY (Z)

Ck(U, Rm)

hình chiếu từ Rn×k lên các không gian

C∞(U, Rm) tập tất cả các ánh xạ trơn

Cω(U, Rm) tập tất cả các ánh xạ giải tích từ U vào Rm

tập tất cả các ánh xạ khả vi liên tục tới cấp k

1

Mở đầu

Với phương trình vi phân (PTVP) ma trận thông thường, nghiệm của chúng là các

hàm ma trận khả vi một biến và thỏa mãn phương trình. Tuy nhiên, khi PTVP đó là

một mô hình toán học mô tả những hiện tượng trong cơ học, hóa học, thiên văn học

thì những ràng buộc đối với ma trận đó thường xuất hiện. Về mặt toán học, chúng đơn

giản chỉ là các biểu thức đại số. Nhưng chúng biểu thị các quy luật bất biến, chẳng

hạn bất biến năng lượng, bất biến về khối lượng, bất biến về thể tích. Trong các ngành

khoa học tương ứng, nó vốn dĩ là các định luật bảo toàn nổi tiếng. Trong rất nhiều

trường hợp, các ràng buộc đó đối với nghiệm đã tạo lên những đa tạp trong không

gian pha. Chẳng hạn, Định lý 1.2.7 trong Chương 1 đã chỉ ra, một đa tạp con nhúng

của Rn có thể được mô tả bởi một biểu thức dạng g(x) = c, trong đó g : Rn → Rm.

Do đó, ta gọi chung những phương trình loại này là PTVP với ràng buộc đa tạp.

Phương pháp giải số các PTVP với ràng buộc đa tạp, ngoài việc đảm bảo các yêu

cầu về tính chính xác và tính ổn định thông thường của một phương pháp số còn phải

đảm bảo rằng nghiệm xấp xỉ luôn nằm trên đa tạp đó. Cụ thể hơn giả sử x0 là điều

xk nói chung ta đều phải có g(xk) = c. Lẽ đương nhiên, các phương pháp được trình

kiện ban đầu của phương trình, tức là g(x0) = c thì với mỗi nghiệm xk tại các bước,

bày ở Chương 1 không đảm bảo việc này. Ta cần phải có những điều chỉnh, thay đổi

thích hợp trong quá trình giải số. Những điều chỉnh này đôi khi tương đối đơn giản

nhưng trong nhiều hoàn cảnh khác lại đòi hỏi những kiến thức sâu sắc về các đa tạp

liên quan.

Với luận văn “Giải số phương trình vi phân ma trận với ràng buộc đa tạp ”, mục

2

tiêu của chúng tôi là trình bày một số cách tiếp cận giải số cho phương trình vi phân

ma trận với ràng buộc đa tạp. Để đạt được điều đó, chúng tôi bố cục luận văn như sau.

Chương 1 được giành để trình bày kiến thức chuẩn bị. Trước tiên, chúng tôi tóm lược

một số phương pháp giải số phương trình vi phân thường phổ biến như phương pháp

BDF, Nystr¨om, Runge-Kutta. Tiếp đó, chúng tôi nhắc lại một số kiến thức về đa tạp

Riemann. Do việc trình bày bài bản tương đối dài, chúng tôi chỉ tóm lược mang tính

diễn giải. Chương 2 là chương chính của luận văn. Trước tiên, chúng tôi chính xác

hóa khái niệm bằng việc trình bày khái niệm bất biến tích phân, lấy ví dụ minh họa.

Sau đó, chúng tôi trình bày phương pháp chính, là phương pháp chung nhất để giải

số phương trình vi phân với ràng buộc đa tạp. Thêm vào đó, chúng tôi cũng trình bày

phương pháp giải số cho hai loại ràng buộc cụ thể trên đa tạp Stiefel và nhóm Lie.

Chúng tôi cũng cố gắng đưa ra một số ví dụ số để minh họa cho các vấn đề được trình

bày. Cuối cùng, luận văn kết thúc bởi phần kết luận và tài liệu tham khảo.

Mặc dù đã rất nghiêm túc và cố gắng thực hiện luận văn này, nhưng luận văn sẽ

khó tránh khỏi những khiếm khuyết nhất định. Kính mong sự góp ý của các thầy cô

để luận văn này được hoàn chỉnh và ý nghĩa hơn.

Luận văn này được thực hiện tại Trường Đại học Khoa học - Đại học Thái Nguyên

và hoàn thành dưới sự hướng dẫn của TS. Nguyễn Thanh Sơn. Tác giả xin được bày

tỏ lòng biết ơn chân thành và sâu sắc tới người hướng dẫn khoa học của mình, người

đã đặt vấn đề nghiên cứu, dành nhiều thời gian hướng dẫn tận tình và đầy trách nhiệm

để tác giả hoàn thành luận văn này.

Tác giả đã học tập được rất nhiều kiến thức chuyên ngành bổ ích cho công tác

và nghiên cứu của bản thân. Nhân dịp này tác giả xin gửi lời cảm ơn sâu sắc tới các

Thầy giáo, Cô giáo đã tham gia giảng dạy lớp Cao học Toán K11C; Nhà trường và các

phòng chức năng của Trường, Khoa Toán - Tin, trường Đại học Khoa học- Đại học

Thái Nguyên đã quan tâm và giúp đỡ tác giả trong suốt thời gian học tập tại trường.

Cuối cùng tác giả xin cảm ơn gia đình, bạn bè, đồng nghiệp đã động viên, ủng hộ

3

và tạo mọi điều kiện cho tác giả trong suốt thời gian nghiên cứu và học tập.

Thái Nguyên, tháng 5 năm 2019

Tác giả luận văn

Đặng Mạnh Hùng

4

Chương 1

Kiến thức chuẩn bị

Trong chương này chúng tôi sẽ trình bày hai mảng kiến thức riêng biệt vốn được

sử dụng trong phần chính của luận văn. Mục đầu tiên của chương được dành để nhắc

lại một vài phương pháp giải số phương trình vi phân quen thuộc. Tài liệu chính cho

mục này là Chương II của quyển sách [5]. Mục thứ hai của chương là một số kiến thức

liên quan đến đa tạp. Để cho tiện, chúng tôi tham khảo hai tài liệu chính [1] và [2].

1.1 Một số phương pháp số giải phương trình vi phân.

Trong mục này, chúng tôi xét phương trình vi phân thường cấp một phụ thuộc vào

˙y = f (t, y), y (t0) = y0.

thời gian

(1.1)

f (t, y1, ..., yn) là một hàm véctơ n thành phần f = (f1, ..., fn). Ta giả sử rằng hệ (1.1)

Ở đây, ta hiểu y = y(t) ∈ Rn là một véctơ cỡ n. Hàm f phụ thuộc vào n + 1 biến

thỏa mãn các điều kiện để nó có nghiệm duy nhất. Các lập luận cũng đúng nếu như y

và f là các ma trận có cấp bất kỳ tương thích với nhau.

5

1.1.1 Phương pháp Runge-Kutta

j=1 aij.

Định nghĩa 1.1.1. Cho bi, aij(i, j = 1, . . . , s) là những số thực và cho ci = (cid:80)s

Một phương pháp Runge-Kutta s-bước được tính bằng cách lấy

s (cid:88)

, i = 1, . . . , s

aijkj

ki = f

t0 + cih, y0 + h

j=1

(cid:32) (cid:33)

s (cid:88)

y1 = y0 + h

biki.

i=1

(1.2)

Các hệ số của phương pháp Runge-Kutta s-bước thường được viết dưới dạng bảng

. . . a1s ...

c1 a11 ...

.

Butcher như sau:

cs as1

. . . ass

. . .

b1

bs

(1.3)

Định nghĩa 1.1.2. Phương pháp Runge-Kutta (hay một phương pháp 1-bước tổng

quát) có p bậc, nếu với mọi bài toán đủ trơn (1.1) sai số địa phương y1 − y (t0 + h)

y1 − y (t0 + h) = O (cid:0)hp+1(cid:1) khi h → 0.

thỏa mãn

Để kiểm tra bậc của phương pháp Runge-Kutta chúng ta phải khai triển chuỗi

Taylor của y (t0 + h) và y1 quanh điểm h = 0. Điều này dẫn đến điều kiện đại số cho

bi = 1

hệ số có các bậc 1, 2, 3 tương ứng như sau:

bici = 1/2

cho bậc 1;

bic2

cho bậc 2; thêm vào đó (1.4)

i = 1/3

thêm vào đó

biaijcj = 1/6

i,j

(cid:88) i (cid:88) i (cid:88) i (cid:88) và cho bậc 3.

Người ta có thể tìm các điều kiện để cho phương pháp có độ chính xác cao hơn. Song,

nó chỉ có ý nghĩa về mặt lý thuyết. Trong thực tế, phương pháp bậc không quá 4 là

tương đối tốt để giải quyết các bài toán nảy sinh.

6

Có thể thấy phương pháp Runge-Kutta có độ khái quát hóa tương đối cao. Chẳng

hạn, phương pháp Euler hiện và phương pháp Euler ẩn là một trường hợp riêng của

phương pháp Runge-Kutta tương ứng với bảng sau

1 1 0 , . (1.5) 1 1

Như vậy, đây là các phương pháp Runge-Kutta một bước. Trong khi đó, phương

pháp hình thang ẩn và phương pháp điểm giữa có thể được mô tả như là một phương

1 suy ra hình thang và trung điểm quy tắc cũng như cả hai phương pháp của Runge

pháp Runge-Kutta bậc hai với các bảng Butcher thứ nhất và thứ hai ở dưới đây. Có bậc

.

0 0 0

1 1/2 1/2 1/2 1/2 1 1/2 1/2 1

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

có bậc 2. Thực tế tính toán cho thấy, phương pháp Kutta bậc 4 dưới đây là thành công

nhất:

0 0

.

1/2 1/2 1/3 1/3

(1.6) 1/2 0 1/2 2/3 -1/3 1

1 0 0 1 -1 1 1 1

1/6 2/6 2/6 1/6 1/8 3/8 3/8 1/8

1.1.2 Phương pháp Runge-Kutta phân hoạch

˙y = f (y, z),

˙z = g(y, z),

Trong mục này chúng tôi xét phương trình vi phân dạng phân hoạch

(1.7)

trong đó y và z có thể là véctơ với kích thước khác nhau. Mục đích là lấy hai phương

pháp Runge-Kutta khác nhau, để xử lý biến y với phương pháp thứ nhất (aij, bi), và

(cid:17) biến z với phương pháp thứ hai . (cid:16) (cid:98)aij,(cid:98)bi

7

Định nghĩa 1.1.3. Cho bi, aij và (cid:98)bi, (cid:98)aij là hệ số của hai phương pháp Runge-Kutta. Phương pháp Runge-Kutta phân hoạch tìm nghiệm của (1.7) được tính bằng cách lấy

s (cid:88)

s (cid:88)

,

y0 + h

aijkj, z0 + h

ki = f

(cid:33) (cid:32)

j=1

j=1

(cid:98)aij(cid:96)j

s (cid:88)

s (cid:88)

,

(cid:33) (cid:32)

y0 + h

aijkj, z0 + h

(cid:96)i = g

j=1

j=1 s (cid:88)

s (cid:88)

y1 = y0 + h

biki, z1 = z0 + h

(1.8) (cid:98)aij(cid:96)j

i=1

i=1

(cid:98)bi(cid:96)i.

Lý thuyết của phương pháp Runge-Kutta có thể được mở rộng một cách trực tiếp

cho phương pháp phân hoạch. Khi (1.8) là phương pháp 1 bước (y1, z1) = Φh (y0, z0),

Định nghĩa 1.1.2 về bậc có thể được áp dụng trực tiếp. Xét bài toán ˙y = f (y), ˙z = g(z),

chúng ta nhìn thầy rằng bậc của (1.8) không thể vượt hơn min(p, (cid:98)p), khi p và (cid:98)p là bậc của hai phương pháp.

Điều kiện cho bậc hai. Khai triển nghiệm chính xác của (1.7) và giải bằng số (1.8)

thành chuỗi Taylor chúng ta thấy rằng phương pháp là bậc 2 nếu điều kiện ghép

bi(cid:98)aij = 1/2,

ij

ij

(cid:88) (cid:88) (1.9) (cid:98)biaij = 1/2,

được thoả mãn bên cạnh các điều kiện cho phương pháp Runge-Kutta bậc hai thông

thường. Chúng tôi cũng lưu ý rằng (1.9) là tự động được thoả mãn bằng phương pháp

phân hoạch do,

ci = (cid:98)ci với mọi y,

(1.10)

j aij và (cid:98)ci = (cid:80)

j (cid:98)aij.

trong đó, ci = (cid:80)

1.1.3 Phương pháp Nystr¨om

¨y = g(t, y, ˙y),

Phương trình vi phân cấp hai

(1.11)

8

tạo thành lớp các bài toán quan trọng. Rất nhiều các bài toán thực tế được mô hình

hóa bởi phương trình vi phân cấp hai (1.11) như bài toán Kepler, hệ thống sử dụng

năng lượng mặt trời ngoài, bài toán trong động lực học phân tử. Lí do chính là nhiều

bài toán cơ học sử dụng Định luật Newton II vốn liên quan đến đạo hàm cấp hai. Đưa

vào biến số mới z = ˙y cho đạo hàm thứ nhất, bài toán (1.11) trở nên tương đương với

˙y = z,

˙z = g(t, y, z).

hệ phân hoạch

(1.12)

s (cid:88)

ki = z0 + h

Một phân hoạch phương pháp Runge-Kutta (1.8) áp dụng cho cho hệ có dạng

j=1

(cid:98)aij(cid:96)j,

s (cid:88)

s (cid:88)

,

(cid:96)i = g

t0 + cih, y0 + h

aijkj, z0 + h

(cid:32) (cid:33)

j=1

j=1

s (cid:88)

s (cid:88)

y1 = y0 + h

biki, z1 = z0 + h

(1.13) (cid:98)aij(cid:96)j

i=1

i=1

(cid:98)bi(cid:96)i.

Nếu chúng ta chèn công thức tính ki vào các phương trình còn lại, ta thu được điều

s (cid:88)

s (cid:88)

aij =

kiện và định nghĩa sau đây, chúng tôi đạt được Định nghĩa 1.1.4 với

aik(cid:98)akj, bi =

bk(cid:98)aki.

k=1

k=1

(1.14)

Định nghĩa 1.1.4. Cho ci, bi, aij và (cid:98)bi, (cid:98)aij là các hệ số thực. Một phương pháp Nystr¨om cho nghiệm của (1.11) được tính bằng cách lấy

s (cid:88)

s (cid:88)

,

(cid:96)i = g

t0 + cih, y0 + cih ˙y0 + h2

aij(cid:96)j, ˙y0 + h

(cid:32) (cid:33)

j=1

j=1

(cid:98)aij(cid:96)j

s (cid:88)

s (cid:88)

y1 = y0 + h ˙y0 + h2

˙y1 = ˙y0 + h

bi(cid:96)i,

(1.15)

i=1

i=1

(cid:98)bi(cid:96)i.

Trong trường hợp đặc biệt quan trọng ¨y = g(t, y), trong đó trường véctơ không

phụ thuộc vào vận tốc, các hệ số aij không cần phải xác định. Phương pháp Nystr¨om là bậc p, nếu y1 − y (t0 + h) = O (cid:0)hp+1(cid:1) và ˙y1 − ˙y (t0 + h) = O (cid:0)hp+1(cid:1).

9

1.2 Khái niệm đa tạp

Trong mục này, chúng tôi sẽ trình bày những khái niệm, định lý, mệnh đề liên

quan đến đa tạp một cách khoa học và cơ bản nhất. Những nội dụng này nhằm bổ trợ

cho các vấn đề ở Chương 2, và được lấy từ tài liệu [2].

1.2.1 Đa tạp khả vi

Định nghĩa 1.2.1. Giả sử (M, τ ) là không gian tô pô Hausdorff với một cơ sở đếm

được. Khi đó M được gọi là một đa tạp tô pô m− chiều nếu nó đồng phôi địa phương

với không gian Rm, nghĩa là với mỗi điểm x ∈ M, tồn tại một lân cận U của x, có một

tập con mở V ⊂ Rm và một phép đồng phôi ϕ : U → V .

Cặp (U, ϕ) được gọi là một bản đồ địa phương hay gọi tắt là bản đồ trong M.

Ta viết Mm để thể hiện đa tạp M có m chiều.

Với U là một tập mở bất kỳ của Rn, ta nhắc lại các kí hiệu Ck(U, Rm), C∞(U, Rm),

và Cω(U, Rm) lần lượt là tập tất cả các ánh xạ khả vi liên tục tới cấp k, tập tất cả các

ánh xạ trơn và tập tất cả các ánh xạ giải tích từ U vào Rm.

Định nghĩa 1.2.2. Xét đa tạp tô pô Mm. Họ A = {(Ui, ϕi) : i ∈ I} các bản đồ trên M

được gọi là một atlas lớp Ck(k ≥ 1) hay Ck− atlas nếu hai điều kiện sau được thỏa

mãn:

(i) họ {Ui} là một phủ mở của M;

i

(ii) với hai bản đồ (Ui, ϕi) và (Uj, ϕj) mà Ui ∩ Uj (cid:54)= ∅ thì ánh xạ chuyển tiếp ϕj ◦ ϕ−1

xác định trên ϕi(Ui ∩ Uj) là ánh xạ khả vi lớp Ck từ ϕi(Ui ∩ Uj) lên ϕj(Ui ∩ Uj).

10

Một bản đồ (U, ϕ) được gọi là tương thích với Ck− atlas A nếu hợp A ∪ {(U, ϕ)} là

một Ck− atlas.

Atlas ˆA được gọi là atlas cực đại nếu nó chứa tất cả các bản đồ tương thích với

nó. Khi đó ˆA cũng được gọi là một Ck− cấu trúc trên Mm.

Cặp (M, ˆA) được gọi là một Ck− đa tạp hay đa tạp khả vi lớp Ck.

Ck− cấu trúc trên M.

Nhận xét 1.2.3. Một Ck− atlas A trên một đa tạp tô pô M xác định duy nhất một

Phát biểu sau đây xây dựng tích Descartes của hai đa tạp.

M1 × M2 là tích Descartes của hai không gian tô pô. Khi đó, tồn tại một atlas A trên

M sao cho (M, ˆA) là một đa tạp khả vi lớp Ck có số chiều

dim M = dim M1 + dim M2.

Mệnh đề 1.2.4. Cho (M1, ˆA1) và (M2, ˆA2) là hai đa tạp khả vi lớp Ck. Cho M =

(M2, ˆA2).

Đa tạp (M, ˆA) được gọi là đa tạp tích Descartes của hai đa tạp (M1, ˆA1) và

1.2.2 Đa tạp con

Định nghĩa 1.2.5. Cho m ≤ n là các số nguyên dương và (N n, ˆAN ) là một đa tạp

khả vi lớp Ck. Một tập con M của N được gọi là một đa tạp con của N nếu với mỗi

11

ϕx(Ux ∩ M ) = ϕx(Ux) ∩ (Rm × {0}),

điểm x ∈ M, tồn tại một bản đồ (Ux, ϕx) ∈ ˆAN sao cho x ∈ Ux và ϕx : Ux ⊂ N → Rm × Rn−m thỏa mãn

số n − m được gọi là đối chiều của M trong N .

Phát biểu sau đây cung cấp chi tiết hơn về cấu trúc khả vi trên đa tạp con.

Mệnh đề 1.2.6. Cho m ≤ n là các số nguyên dương và (N n, ˆAN ) là một đa tạp khả

π : Rm × Rn−m → Rm

vi lớp Ck. Cho M là một đa tạp con của N và được trang bị tô pô con. Ký hiệu

là phép chiếu lên thành phần thứ nhất. Khi đó

AM := {(Ux ∩ M, (π ◦ ϕx) |Ux∩M) : x ∈ M}

là một atlas lớp Ck trên M. Vì thế (M, ˆAM) là một đa tạp khả vi m chiều lớp Ck.

Cấu trúc khả vi ˆAM trên M xác định trong mệnh đề trên được gọi là cấu trúc cảm

sinh từ ˆAN .

min{m, n}, một ánh xạ tuyến tính giữa các không gian hữu hạn chiều đều có thể coi

Ta nhắc lại rằng một ma trận A cỡ m × n được gọi là hạng đủ nếu rankA =

tương đương như một ma trận. Định lý sau đây cho ta một nguồn dồi dào những ví dụ

về đa tạp con.

df |x : Rn → Rm

Định lí 1.2.7 (Ánh xạ ẩn). Cho m ≤ n là các số nguyên dương và f : U → Rm là một ánh xạ lớp Ck từ một tập mở U ⊂ Rn. Nếu x ∈ U, f (x) = y và Jacobian của nó

là hạng đủ thì f −1({y}) là một đa tạp con khả vi lớp Ck của Rn có số chiều n − m.

12

1.2.3 Vectơ tiếp xúc, không gian tiếp xúc

Trước tiên, ta định nghĩa khái niệm này trong không gian Rn quen thuộc.

Định nghĩa 1.2.8. Cho x là một điểm trong Rm và ký hiệu TxRm là tập các toán tử vi phân tuyến tính tại x triệt tiêu hằng số. Tức là TxRm gồm các ánh xạ α : ε(x) → R

thỏa mãn

(i) α(λf + µg) = λα(f ) + µα(g),

(ii) α(f g) = α(f ) · g(x) + f (x)α(g),

với mọi α, µ ∈ R và f, g ∈ ε(x).

Dễ dàng nhận thấy tập TxRm có cấu trúc của một không gian véctơ thực với hai

(α + µ)(f ) := α(f ) + β(f ),

(λα)(f ) := λα(f ).

phép toán

φ : Rm → TxRm

v (cid:55)→ ∂v

Định lí 1.2.9. Cho x ∈ Rm. Khi đó, ánh xạ

là một đẳng cấu giữa các không gian véctơ.

Bây giờ ta xét các khái niệm đó đối với đa tạp.

Định nghĩa 1.2.10. Cho M là một đa tạp khả vi, x ∈ M và ε(x) là tập các hàm

ξx : ε(x) → R thỏa mãn

thực định nghĩa trên một lân cận mở của x. Một véctơ tiếp xúc ξx tại x là một ánh xạ

(i) ξx(λf + µg) = λξx(f ) + µξx(g),

13

(ii) ξx(f g) = ξx(f )g(x) + f (x)ξx(g),

với mọi λ, µ ∈ R và f, g ∈ ε(x).

Tập các vectơ tiếp xúc của M tại x được gọi là không gian tiếp xúc tại x và ký

hiệu TxM.

(ξx + ζx)(f ) = ξx(f ) + ζx(f ),

(λξx)(f ) = λξx(f ),

Phép cộng và phép nhân với vô hướng trong TxM được định nghĩa như sau

với mọi ξx, ζx ∈ TxM, f ∈ ε(x) và λ ∈ R.

1.3 Đa tạp Riemann

1.3.1 Khái niệm

Cho M là một đa tạp trơn, tập các hàm trơn C∞(M) trên M lập thành một vành

giao hoán. Tập các hàm trơn C∞(T M) trên M lập thành một môđun trên C∞(M).

0 (T M) = C∞(M) và với số nguyên dương r, đặt

C∞ r (T M) = C∞(T M) ⊗ ... ⊗ C∞(T M),

Đặt C∞

trong đó ⊗ là ký hiệu tích tenxơ của các không gian vectơ.

Định nghĩa 1.3.1. Cho M là một đa tạp khả vi. Một mêtric Riemann g trên M là một

2 (T M) → C∞

0 (T M) sao cho với mỗi x ∈ M, hạn chế gx

ánh xạ song tuyến tính : C∞

gx(ξx, ζx) = g(ξ, ζ)x = (cid:104)ξx, ζx(cid:105)x.

của g trên TxM ⊗ TxM với

là một tích vô hướng trên không gian tiếp xúc TxM. Cặp (M, g) được gọi là một đa

tạp Riemann.

14

Định nghĩa 1.3.2. Cho M là một đa tạp con của đa tạp Riemann (M, g). Khi đó mọi

không gian tiếp xúc TxM đều có thể coi như là một không gian con của TxM. Mêtric

ξ, ζ, ∈ TxM,

gx(ξ, ζ) := gx(ξ, ζ),

Riemann trong M có thể được định nghĩa như sau

trong đó ξ và ζ ở vế phải được xem như là các phần tử của TxM. Được trang bị mêtric

như trên, M trở thành một đa tạp Riemann, gọi là đa tạp Riemann con của M.

M tại x:

(TxM)⊥ = {ξ ∈ TxM : gx(ξ, ζ) = 0, ∀ ζ ∈ TxM}.

Phần bù trực giao của TxM trong TxM được gọi là không gian pháp tuyến của

Mọi phần tử ξ ∈ TxM đều có thể biểu diễn duy nhất dưới dạng tổng của một phần tử

ξ = Pxξ + P ⊥

x ξ,

của TxM và một phần tử của (TxM)⊥:

x là phép chiếu trực giao lên

(TxM)⊥.

trong đó Px là phép chiếu trực giao lên TxM và P ⊥

Mêtric Riemann cho phép ta lượng hóa, đo đạc trên các không gian tiếp xúc. Hơn

thế nữa, nó giúp xây dựng cấu trúc không gian mêtric từ cấu trúc hình học trên đa tạp

Riemann như khái niệm khoảng cách trong mục sau đây.

1.3.2 Khoảng cách

C1 trong M. Khi đó, độ dài của γ được định nghĩa là

Định nghĩa 1.3.3. Cho (M, g) là một đa tạp Riemann và γ : I → M là một cung lớp

L(γ) =

I

(cid:90) (cid:112)g( ˙γ(t), ˙γ(t))dt. (1.16)

x, y ∈ M, ký hiệu CM là tập tất cả các cung γ : [0, 1] → M sao cho γ(0) = x, γ(1) = y

Mệnh đề 1.3.4. Cho (M, g) là một đa tạp Riemann liên thông đường. Với mỗi điểm

15

d : M × M → R+

(x, y)

(cid:55)→ inf{L(γ) : γ ∈ CM}.

và định nghĩa hàm

Khi đó, (M, d) lập thành một không gian mêtric. Hơn nữa, tô pô sinh bởi d trùng với

tô pô ban đầu của M.

Ta có định lý quan trọng sau đây về sự tồn tại của mêtric Riemann.

Định lí 1.3.5. Cho (Mm, ˆA) là một đa tạp khả vi. Khi đó, tồn tại một mêtric Riemann

trên M.

16

Chương 2

Bảo toàn của tích phân và phương pháp giải

Trong chương này, chúng tôi trình bày một cách chi tiết cách giải các phương trình

vi phân ma trận với ràng buộc đa tạp. Những nội dung này chủ yếu dựa vào tài liệu [3]

và Chương IV tài liệu [5].

2.1 Khái niệm bất biến trong tích phân

˙y = f (y),

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

(2.1)

trong đó y có thể là một véctơ hoặc một ma trận. Thực ra, một véctơ cũng có thể coi

là một ma trận có một cột nếu ta có thể coi (2.1) một cách tổng quát là một phương

trình vi phân đối với ma trận.

Định nghĩa 2.1.1. Một hàm khác hằng số I(y) được gọi là một tích phân đầu (first

I (cid:48)(y)f (y) = 0, với mọi y.

integral) của phương trình (2.1) nếu

(2.2)

I (y(t))(cid:48) = I (cid:48)(y) ˙y = I (cid:48)(y)f (y) = 0.

Từ (2.2), ta nhận thấy ngay

Nói cách khác, giá trị của hàm I dọc theo nghiệm y(t) của phương trình (2.1) là hằng

số. Do tính chất này mà tích phân đầu còn được gọi là một bất biến, một đại lượng

17

bảo toàn (conserred quantity) hay một chuyển động hằng (constant motion). Ngay

sau đây, ta sẽ xét một ví dụ để mình họa khái niệm này.

Ví dụ 2.1.2. (Hệ Hamilton)

Cho H = H(p, q) = H(p1, . . . , pn, q1, . . . , qn) là hàm 2n biến số phụ thuộc thời gian t

˙p = −Hq(p, q), ˙q = Hp(p, q),

(phương trình vi phân) Hamilton là hệ thỏa mãn

(2.3)

trong đó Hp = ∇pH, Hq = ∇qH, tức là các gradiant theo từng biến số. Xét đạo hàm

H (cid:48) (p(t), q(t)) = ∇pH T ˙p + ∇qH T ˙q

= ˙q ˙p − ˙p ˙q

= 0.

của H theo t

Do đó, H là hàm hằng theo t. Trong trường hợp p(t) là động lượng, q(t) = x(t) là tọa

p2 +

x2,

H(p, x) =

độ của vật thì hàm

1 2

1 2

(2.4)

thể hiện năng lượng toàn phần của hệ. Có thể kiểm tra thấy H ở (2.4) thỏa mãn hệ

(2.3). Do H là hằng theo t nên ta nói năng lượng toàn phần của hệ được bảo toàn. Khi

giải số phương trình vi phân (2.3), đại lượng (2.4) phải được bảo toàn tại mọi bước

thời gian.

Ví dụ 2.1.3. (Bảo toàn khối lượng trong phản ứng hóa học)

Ta xét bài toán Robertson sau đây. Giả sử có 3 chất A, B, C tham gia một phản ứng

0.04 −−→ B,

A

B + B

3×107 −−−−→ C + B,

B + C

104 −−→ A + C.

hóa học với tốc độ hằng như sau

18

B, C lần lượt là y1, y2, y3. Nguyên lý phản ứng hóa học phát biểu rằng tốc độ biến

Ký hiệu nồng độ (do thể tích không đổi nên đồng thời cũng là lượng) của các chất A,

thiên của chất sẽ tỉ lệ với chất tạo thành bởi hằng số (gọi là tốc độ phản ứng). Theo

đó, ta có hệ phương trình

˙y1 = −0.04y1 + 104y2y3



˙y2 = 0.04y1 − 104y2y3 − 3 · 107y2 2

(2.5)

˙y3 = 3 · 107y2 2.

 

˙y1 + ˙y2 + ˙y3 = 0,

Có thể thấy ngay nghiệm của hệ phương trình này thỏa mãn

y1 + y2 + y3 = const.

hay

Điều này trên thực tế phản ánh Định luật bảo toàn khối lượng trong phản ứng hóa học.

Đây là một đại lượng bất biến cần được bảo toàn khi ta tích phân số hệ phương trình

vi phân (2.5).

Sau đây, ta sẽ xét kiểu bất biến tích phân đơn giản.

I(y) = d(cid:62)y, trong đó d là một véctơ hoặc ma trận hằng.

Định nghĩa 2.1.4. Bất biến I(y) được gọi là một bất biến tuyến tính nếu nó có dạng

Định lí 2.1.5. (Bảo toàn bất biến tuyến tính) Tất cả phương pháp Runge-Kutta ẩn và

hiện đều bảo toàn bất biến tuyến tính. Phương pháp Runge-Kutta phân hoạch bảo

toàn bất biến tuyến tính nếu bi = (cid:98)bi với mọi i, hoặc nếu bất biến chỉ phụ thuộc vào p

hoặc q.

Chứng minh. Đặt I(y) = d(cid:62)y với một véctơ hằng d, sao cho d(cid:62)f (y) = 0 với mọi y.

Trong trường hợp của phương pháp Runge-Kutta chúng ta luôn có d(cid:62)ki = 0, và vì vậy

19

d(cid:62)y1 = d(cid:62)y0 + hd(cid:62) (cid:0)(cid:80)s

i=1 biki

(cid:1) = d(cid:62)y0. Chứng minh cho phương pháp phân hoạch

cũng tương tự như vậy.

˙y = A(y)y,

Để kết thúc mục này, ta xét phương trình vi phân có dạng đặc biệt

(2.6)

trong đó y có thể là ma trận không nhất thiết vuông.

Định lí 2.1.6. Nếu A(y) là ma trận phản đối xứng, tức là A(y)(cid:62) = −A(y) thì hàm bậc

hai I(y) = y(cid:62)y là một bất biến. Đặc biệt, nếu giá trị ban đầu y0 có các cột trực chuẩn

thì các cột của nghiệm y(t) của (2.6) cũng luôn trực chuẩn với mọi t.

I (cid:48)(Y )f (Y ) = I (cid:48)(Y ) (A(Y )Y ) = Y (cid:62)A(Y )Y + Y (cid:62)A(Y )(cid:62)Y triệt tiêu với mọi Y , vì A(Y )

Chứng minh. Đạo hàm của I(y) là I (cid:48)(Y )H = Y (cid:62)H + H (cid:62)Y . Vì vậy, chúng ta có

là phản đối xứng. Phát biểu được chứng minh.

Ví dụ 2.1.7. (Vật rắn không biến dạng (rigid body))

Rigid body là thuật ngữ để chỉ một vật mà khoảng cách giữa hai điểm bất kỳ của

vật là không thay đổi. Người ta đã chỉ ra chuyển động của vật rắn không biến dạng tự

˙y1 = a1y2y3, a1 = (I2 − I3) / (I2I3)

do thỏa mãn phương trình Euler:

˙y2 = a2y3y1, a2 = (I3 − I1) / (I3I1)

˙y3 = a3y1y2, a3 = (I1 − I2) / (I1I2) ,

(2.7)

trong đó, y = (y1, y2, y3)(cid:62) là mômen động lượng (argular momentum) và I1, I2, I3 là

các mômen quán tính. Từ (2.7), ta có thể viết lại hệ phương trình dưới dạng

0

˙y1

y3/I3 −y2/I2

y1

=

.

     

0

y2

˙y2

−y3/I3

y1/I1

(2.8)

0

˙y3

y2/I2 −y1/I1

y3

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

20

Đến đây, phương trình (2.8) có dạng (2.6) với ma trận hệ số thỏa mãn điều kiện phản

1 + y2

2 + y2

3 là một bất biến. Trên thực

đối xứng của Định lý 2.1.6. Theo định lý này, y2

+

+

,

H(y1, y2, y3) =

1 2

tế (cid:19)

y2 2 I2

y2 3 I3

(cid:18)y2 1 I1

là một bất biến bậc hai mô tả động năng của hệ.

2.2 Bất biến bậc hai

Q(y) = y(cid:62)Cy,

Trong mục này, ta sẽ xét bất biến bậc hai dạng

trong đó, C là một ma trận đối xứng. Cụ thể, nó là một bất biến nếu y(cid:62)Cf (y) = 0. Ta

sẽ lần lượt khảo sát xem dưới điều kiện nào thì phương pháp Runge-Kutta và phương

pháp Nystr¨om bảo toàn các bất biến bậc hai.

biaij + bjaji = bibj, ∀i, j = 1, · · · , s,

Định lí 2.2.1. Nếu hệ số của phương pháp Runge-Kutta thoả mãn

(2.9)

thì nó là bảo toàn bất biến bậc hai.

i=1 biki của Định nghĩa 1.1.3 cho

s (cid:88)

s (cid:88)

s (cid:88)

bik(cid:62)

bjy(cid:62)

bibjk(cid:62)

Chứng minh. Biểu thức y1 = y0 + h (cid:80)s

1 Cy1 = y(cid:62) y(cid:62)

0 Cy0 + h

0 Ckj + h2

i Cy0 + h

i Ckj.

j=1

i=1

i,j=1

(2.10)

j=1 aijkj. Ý tưởng chính là rút y0 từ quan hệ

Ta biểu thị ki = f (Yi) với Yi = y0 + h (cid:80)s

này và đưa nó vào biểu thức trung tâm của (2.10). Sử dụng phép đối xứng của C cho

s (cid:88)

s (cid:88)

biY (cid:62)

(bibj − bjaij − bjaji) k(cid:62)

1 Cy1 = y(cid:62) y(cid:62)

0 Cy0 + h

i Cf (Yi) + h2

i Ckj.

i,j=1

i=1

ta

Điều kiện (2.9) cùng với giả thiết y(cid:62)Cf (y) = 0, mà thực chất phát biểu rằng y(cid:62)Cy là

1 Cy1 = y(cid:62)

0 Cy0.

bất biến của (2.1), suy ra rằng y(cid:62)

21

Sau đây, ta sẽ xét một bất biến bậc hai tổng quát hơn dạng đã xét ở trên. Xét hệ

˙y = f (y, z), ˙z = g(y, z).

phương trình vi phân (hoặc có thể hiểu là phương trình vi phân đã được phân hoạch):

(2.11)

Q(y, z) = y(cid:62)Dz,

Đối với hệ dạng (2.11), xét bất biến bậc hai có dạng

(2.12)

trong đó D là một ma trận hằng. Định lý sau đây có thể được coi là một mở rộng của

Định lý 2.2.1.

Định lí 2.2.2. Nếu hệ số của phương pháp Runge-Kutta phân hoạch thỏa mãn

bi(cid:98)aij + (cid:98)bjaji = bi(cid:98)bj; i, j = 1, · · · , s,

(2.13)

bi = (cid:98)bj; i = 1, · · · , s,

(2.14)

thì nó bảo toàn bất biến bậc hai của dạng (2.12). Nếu phương trình vi phân phân

hoạch đặc biệt có dạng ˙y = f (z), ˙z = g(y) thì chỉ cần điều kiện (2.13) để suy ra rằng

bất biến dạng (2.12) được bảo toàn.

s (cid:88)

s (cid:88)

s (cid:88)

bik(cid:62)

y(cid:62) 1 Dz1 = y(cid:62)

0 Dz0 + h

0 D(cid:96)j + h2

i Dz0 + h

i D(cid:96)j.

Chứng minh. Chứng minh gần giống với Định lý 2.2.1 thay vì (2.10), ta có

bi(cid:98)bjk(cid:62)

i=1

j=1

i,j=1

(cid:98)bjy(cid:62)

Gọi (Yi, Zi) là các biến của hàm ki = f (Yi, Zi) và (cid:96)j = g (Yi, Zi), sử dụng kỹ thuật

s (cid:88)

s (cid:88)

bif (Yi, Zi)(cid:62) DZi + h

1 Dz1 = y(cid:62) y(cid:62)

0 Dz0 + h

j Dg (Yj, Zj)

giống như trong Định lý 2.9 cho

j=1

(cid:98)bjY (cid:62)

i=1 s (cid:88)

+ h2

k(cid:62) i D(cid:96)j.

bi(cid:98)bj − bi(cid:98)aij − (cid:98)bjaji

i,j=1

(cid:16) (cid:17) (2.15)

Vì (2.12) là một bất biến, chúng ta có f (y, z)(cid:62)Dz + y(cid:62)Dg(y, z) = 0 với mọi y và z. Vì

1 Dz1 = y(cid:62)

0 Dz0.

vậy hai điều kiện (2.13) và (2.14) dẫn đến y(cid:62)

22

Với trường hợp đặc biệt trong đó f chỉ phụ thuộc vào z và g chỉ phụ thuộc vào y,

f (z)(cid:62)Dz = −y(cid:62)Dg(y) = const.

giả thiết f (z)(cid:62)Dz + y(cid:62)Dg(y) = 0 (cho tất cả y, z) suy ra

Do đó, điều kiện (2.14) không còn cần thiết cho việc chứng minh nữa.

¨y = g(y),

Có một điểm thú vị là phương trình vi phân cấp hai dạng

(2.16)

˙y = z, ˙z = g(y).

có thể được xem là dạng đặc biệt của (2.11).

Như đã biết ở Mục 1.1.3, phương pháp Nystr¨om được thiết kế chuyên để giải phương

trình loại này. Ta hãy tìm hiểu xem với điều kiện nào thì phương pháp này bảo toàn

bất biến dạng (2.12).

βi =

bi (1 − ci) với i = 1, · · · , s,

Định lí 2.2.3. Nếu hệ số của phương pháp Nystr¨om thỏa mãn

bi (βj − aij) = bj (βi − aji) với i, j = 1, · · · , s,

(2.17)

thì nó bảo toàn bất biến bậc hai dạng y(cid:62)D ˙y.

Chứng minh. Dạng toàn phương Q(y, ˙y) = y(cid:62)D ˙y là một tích phân đầu của (2.16) nếu

˙y(cid:62)D ˙y + y(cid:62)Dg(y) = 0, với mọi y, ˙y ∈ Rn.

và chỉ nếu

(2.18)

Điều này suy ra D là đối xứng lệch và y(cid:62)Dg(y) = 0. Giống như chứng minh của Định

1 D ˙y1 sử dụng công thức

lý 2.2.1 và Định lý 2.2.2 bây giờ chúng ta tính y(cid:62)

s (cid:88)

(cid:96)i = g

y0 + cih ˙y0 + h2

aij(cid:96)j

j=1

(cid:32) (cid:33)

s (cid:88)

s (cid:88)

y1 = y0 + h ˙y0 + h2

˙y1 = ˙y0 + h

βi(cid:96)i,

bi(cid:96)i

i=1

i=1

(2.19)

23

j aijlj, trong đó Yi ký hiệu là biến của g ở

và chúng ta thế y0 bởi Yi − cih ˙y0 − h2 (cid:80)

s (cid:88)

biY (cid:62)

1 D ˙y1 = y(cid:62) y(cid:62)

0 D ˙y0 + h ˙y(cid:62)

0 D ˙y0 + h

i D(cid:96)i

i=1 s (cid:88)

s (cid:88)

+ h2

bi (1 − ci) ˙y(cid:62)

βi(cid:96)(cid:62)

0 D(cid:96)i

i D ˙y0 + h2

i=1

i=1 s (cid:88)

+ h3

bi (βj − aij) (cid:96)(cid:62)

j D(cid:96)i.

i,j=1

(2.19). Điều này dẫn đến

i D(cid:96)i = Y (cid:62)

i Dg (Yi) = 0, điều kiện (2.17) suy ra

Dùng tính đối xứng lệch của D và Y (cid:62)

1 D ˙y1 = y(cid:62)

0 D ˙y0.

tính chất bảo toàn y(cid:62)

2.3 Phương pháp chiếu

Trong mục này, ta sẽ xét ràng buộc đa tạp dưới dạng một biểu thức. Như đã biết

ở Định lý 1.2.7 một đa tạp con nhúng n − m chiều của Rn có thể được biểu thị như là

M = {y : g(y) = 0},

một tập thỏa mãn ràng buộc

(2.20)

˙y = f (y),

trong đó g : Rn → Rm có đạo hàm là một toàn ánh. Xét phương trình vi phân

y0 ∈ M ⇒ y(t) ∈ M, ∀y.

với tính chất

(2.21)

Lưu ý rằng điều kiện (2.21) yếu hơn điều kiện (2.2) trong Định nghĩa 2.1.1. Thật vậy,

điều kiện (2.21) có thể được viết lại là g(cid:48)(y)f (y) = 0, với mọi y ∈ M trong khi điều

kiện (2.2) là g(cid:48)(y)f (y) = 0, với mọi y ∈ Rn. Do đó g(y) thỏa mãn (2.21) được gọi là

một bất biến yếu và ta nói rằng phương trình ˙y = f (y) là một phương trình vi phân

trên đa tạp M. Trong trường hợp M là đa tạp con nhúng của Rn, sự tồn tại nghiệm

của phương trình vi phân trên đa tạp M được đảm bảo bởi định lý sau đây.

24

˙y = f (y) là một phương trình vi phân trên đa tạp M nếu và chỉ nếu

f (y) ∈ M với mọi y ∈ M.

Định lí 2.3.1. Giả sử M là một đa tạp con của Rn. Khi đó, phương trình vi phân

Ví dụ 2.3.2. Xét phương trình con lắc đơn trong không gian 2 chiều với tọa độ

˙q1 = p1, ˙p1 = −q1λ

Descartes.

˙q2 = p2, ˙p2 = −1 − q2λ,

(2.22)

1 + p2

2 − q2)/(q2

1 + q2

2). Xét

H1(p, q) = p(cid:62)q = p1q1 + p2q2.

với λ = (p2

=

∂H1(p, q) ∂t

∂(p1q1 + p2q2) ∂t

= p1

+ q1

+ p2

+ q2

dq1 dt

dp1 dt

dq2 dt

dp2 dt

= p1 ˙q1 + q1 ˙p1 + p2 ˙q2 + q2 ˙p2.

Ta có

= p1p1 + q1(−q1λ) + p2p2 + q2(−1 − q2λ)

∂H1(p, q) ∂t

= p2

1 + p2

2 − q2 − λ(q2

1 + q2 2).

Theo Giả thiết 2.22 nên

1 + p2

2 − q2)/(q2

1 + q2

2) nên suy ra

= p2

1 + p2

2 − q2)/(q2

1 + q2

2)(cid:3) (q2

1 + q2 2)

1 + p2

2 − q2 − (cid:2)(p2

∂H1(p, q) ∂t

= p2

2 + q2

1 − p2

2 − q2 − p2

1 + p2

= 0.

Mặt khác λ = (p2

= 0. Như vậy, đây là một bất biến (theo Định nghĩa 2.1.1 được thảo

∂H1(p, q) ∂t luận ở phần trên).

Vậy

25

1 + q2

2. Ta có

=

(q2

∂H2(p, q) ∂t

∂ ∂t

+ 2q2

= 2q1

1 + q2 2) dq1 dt

dq2 dt

= 2(q1p1 + p2q2) = 2H1(p, q).

Tuy nhiên xét H2(p, q) = q2

= 0 nên

∂H1(p, q) ∂t

= 0.

H2(p, q) = 2

∂H1(p, q) ∂t

Mà H1(p, q) không phụ thuộc thời gian t tức là

Điều đó nghĩa là H2 là một bất biến yếu.

Khi giải số những phương trình loại này bằng một trong các phương pháp trình

bày ở Chương 1, nhìn chung các bất biến kể trên không được đảm bảo. Cụ thể hơn, với

giả thiết y0 ∈ M, ta tính ˜y1 bởi một phương pháp nào đó theo y0, ký hiệu ˜y1 = φh(y0)

trong đó h là độ dài bước thời gian. Nhìn chung ˜y1 /∈ M. Một ý tưởng tự nhiên là thay

vì sử dụng ˜y1, ta sử dụng hình chiếu của nó xuống M như là một xấp xỉ.

y1 : = PM(˜y1) = PM (φh(y0))

(2.23)

Nếu ta đảm bảo được (2.23) trong mọi bước giải, ta sẽ thu được nghiệm xấp xỉ của

phương trình mà ràng buộc đa tạp g(yk) = 0 được thỏa mãn với mọi k.

Phép chiếu PM phụ thuộc vào đa tạp cụ thể khi tính toán. Tuy nhiên nó phải đảm

bảo rằng sai số giữa yk và ˜yk không quá lớn, chẳng hạn O(hp+1) thì phương pháp chiếu

sẽ không ảnh hưởng đến phương pháp giải số Oh có độ chính xác O(hl) với l ≤ p.

Trước khi chuyển sang mục sau, chúng ta thảo luận về tính tồn tại nghiệm của

phương trình vi phân thỏa mãn các bất biến. Trên thực tế, các phương trình vi phân

đã xét là mô hình toán học cho các hiện tượng trong Vật lý, Hóa học, Cơ học như ở

các Ví dụ 2.1.2, Ví dụ 2.1.3, Ví dụ 2.1.7. Những bất biến đó cũng là mô tả dưới dạng

toán học của các quy luật về động lượng, về bảo toàn vật chất hay về cấu trúc cụ thể

của hệ cơ học. Do đó mà sự tồn tại nghiệm của những phương trình vi phân được đề

26

cập (giả sử các mô hình đó là đúng đắn) được xác định bằng thực nghiệm. Trên thực

tế, thảo luận về tính duy nhất nghiệm mới là bài toán thú vị hơn. Qua tính toán, xem

Mục 2.6, bản thân các phương trình vi phân có thể có nhiều hơn một nghiệm do các

ràng buộc vi phân chưa đủ để dẫn tới tính duy nhất. Chính các ràng buộc đại số (các

bất biến) đã giúp cho phương trình vi phân có nghiệm duy nhất nghiệm.

2.4 Tích phân phương trình vi phân trên đa tạp Stiefel

Trong mục này, chúng tôi sẽ trình bày phương pháp tích phân số bảo toàn ràng

buộc trên đa tạp Stiefel. Phần này được viết dựa trên tài liệu [3].

2.4.1 Sơ lược về đa tạp Stiefel

Nhắc lại rằng đa tạp Stiefel V(k, n) là tập các ma trận cỡ n × k (n ≥ k) có các cột

V(k, n) = {Y ∈ Rn×k : Y T Y = I}.

trực chuẩn:

Không gian tiếp xúc của V(k, n) lại Y được xác định bởi

TY V(k, n) = {∆ ∈ Rn×k, Y (cid:62)∆ + ∆(cid:62)Y = 0}.

(2.24)

Với tích vô hướng Euclide của không gian Rn×k chứa V(k, n) như một đa tạp con, thì

NY = {N ∈ Rn×k : (cid:104)∆, N (cid:105) = 0, ∀∆ ∈ TyV(n, k)}.

ta hoàn toàn xác định được không gian chuẩn tắc của V(k, n) lại Y

sym(A) = (cid:0)A + A(cid:62)(cid:1) /2, skew(A) = (cid:0)A − A(cid:62)(cid:1) /2,

Ký hiệu

là phần đối xứng và phản đối xứng của một ma trận vuông A bất kỳ. Hình chiếu từ

Rn×k lên các không gian tiếp xúc và không gian chuẩn tắc của đa tạp được cho bởi

ΠTY (Z) = Yskew(Y (cid:62)Z) + (cid:0)I − Y Y (cid:62)(cid:1) Z,

công thức

27

ΠNY (Z) = Ysym(Y (cid:62)Z).

và

2.4.2 Cung trắc địa trên đạ tạp Stiefel

Mêtric hiện tại được sử dụng trên đa tạp Stiefel V(k, n) là mêtric chính tắc trong

không gian Euclide Rn×k vốn chứa V(k, n) như một đa tạp con nhúng. Sử dụng mêtric

này để định nghĩa khoảng cách trên đa tạp V(k, n) sẽ không phù hợp vì nó vốn dĩ là

khoảng cách thẳng thuần túy trong không gian, trong khi V(k, n) là một đa tạp không

tầm thường. Do đó, ta cần chọn mêtric khác, gọi là mêtric chính tắc và định nghĩa nó

như sau

Y Y (cid:62)(cid:17)

Φ

∆(cid:62) (cid:16)

I −

, ∀∆, Φ ∈ TY V(k, n)

(cid:104)∆, Φ(cid:105)Y = trace

1 2

(cid:17) (cid:16) (2.25)

Tiếp theo, ta xây dựng cung trắc địa tương ứng với mêtric (2.25). Một cung trắc

địa G(t) được hiểu là mở rộng của khái niệm đường thẳng trong không gian lên đa

tạp. Nó là đường liên tục có độ dài ngắn nhất nối 2 điểm với nhau. Cung trắc địa được

xác định khi biết 2 đầu mút hoặc 1 đầu mút cùng với hướng đi ban đầu của nó. Trong

trường hợp thứ nhất việc xác định cung trắc địa tương đương việc giải một bài toán

biên của một phương trình vi phân cấp hai, trong khi trường hợp thứ hai là một bài

toán giá trị ban đầu. Người ta đã chỉ ra (xem [4]) phương trình vi phân xác định cung

trắc địa trên V(k, n) ứng với mêtric (2.25) có dạng

= 0.

+ ˙G(cid:62) ˙G

¨G + ˙G ˙G(cid:62)G + G

(cid:17) (cid:16)(cid:0)G(cid:62) ˙G(cid:1)2 (2.26)

Có thể thấy (2.26) là một phương trình vi phân ma trận phi tuyến nên rất khó có thể

giải được trực tiếp. Ý tưởng chính ở đây là xấp xỉ (2.26) bởi một phương trình đơn

giản hơn. Định lý sau đây cho ta một phát biểu tường minh.

Định lí 2.4.1. Cho Y0 ∈ V(k, n)(R), H ∈ TY0 và xét hàm ma trận sau:

G(t) = Y0M (t) + (cid:0)In − Y0Y (cid:62)

0

(cid:1) HN (t), t > 0, (2.27)

28

˙M (t) = BM (t) + CN (t), M (0) = Ip,

khi M (t), N (t) ∈ Rp×p là giải được hệ phương trình vi phân tuyến tính

˙N (t) = M (t), N (0) = 0,

˙G(0) = H, và G(t)

(2.28)

0

0 H và C = −H (cid:62) (cid:0)In − Y0Y (cid:62)

với B = Y (cid:62) (cid:1) H. Khi đó G(0) = Y0,

luôn nằm trên V(k, n)(R) với bất kì t > 0. Hơn nữa, nếu Y (t) là nghiệm của hệ phương

˙Y (t) = F (Y (t))Y (t), t > 0, Y (0) = Y0 ∈ Vn,p(R),

trình vi phân

Y (t), nghĩa là G(t) − Y (t) = 0(t2).

và chúng ta giả sử rằng Y0 = Y (0) và H = ˙Y (0), thì G(t) là một xấp xỉ thứ nhất của

Chứng minh. Bằng cách xây dựng rằng G(0) = Y0, và vì

˙G(t) = Y0 ˙M (t) + (cid:0)In − Y0Y (cid:62)

0

(cid:1) HM (t), t > 0, (2.29)

G(cid:62)(t)G(t) = ˙G(cid:62)(t)G(t) + G(cid:62)(t) ˙G(t) = 0, t > 0,

d dt

chúng ta đạt được ˙G(0) = H. Bây giờ, G(t) thuộc vào V(k, n)(R) nếu và chỉ nếu

˙G(cid:62)(t)G(t) là đối xứng lệch. Chúng ta chú ý rằng

˙G(cid:62)(t)G(t) = (cid:2) ˙M (cid:62)(t)Y (cid:62)

tức là khi

0

0

(cid:1) HN (t)(cid:3)

0

(cid:1)(cid:3) (cid:2)Y0M (t) + (cid:0)In − Y0Y (cid:62) (cid:1) HN (t)

0 + M (cid:62)(t)H (cid:62) (cid:0)In − Y0Y (cid:62) = ˙M (cid:62)(t)M (t) + M (cid:62)(t)H (cid:62) (cid:0)In − Y0Y (cid:62) = M (cid:62)(t)B(cid:62)M (t) + N (cid:62)(t)C(cid:62)M (t) + M (cid:62)(t)H (cid:62) (cid:0)In − Y0Y (cid:62)

0

(cid:1) HN (t).

0

˙G(cid:62)(t)G(t) = M (cid:62)(t)B(cid:62)M (t) + N (cid:62)(t)CM (t) − M (cid:62)(t)CN (t),

Do B(cid:62) = −B và ma trận C = −H (cid:62) (cid:0)I − Y0Y (cid:62) (cid:1) H là đối xứng, thì cho nên

đó là ˙G(cid:62)(t)G(t) là đối xứng lệnh với mọi t > 0. Bây giờ, sử dụng khai triển Taylor của

29

M (t) và N (t) và điều kiện ban đầu, từ 2.27 ta suy ra

M (0) + ˙t ˙M (0) +

G(t) =Y0

t2 2

(cid:20) (cid:21) ¨M (0)

+ 0 (cid:0)t3(cid:1)

+ (cid:0)In − Y0Y T

0

t2 2

(cid:21) ¨N (0) (cid:1) H (cid:20) N (0) + t ˙N (0) +

Y0

=Y0 + tY0B +

0

t2 2

+

(2.30) (cid:1) H (cid:0)B2 + C(cid:1) + t (cid:0)In − Y0Y T

0

t2 2

=Y0 + tH +

(Y0C + HB) + 0 (cid:0)t3(cid:1)

(cid:1) HB + 0 (cid:0)t3(cid:1)

(cid:0)In − Y0Y T t2 2

Vì vậy, nếu H = ˙Y (0) thì G(t) − Y (t) = 0(t2).

Như vậy, thay cho việc giải phương trình vi phân ma trận (2.26), ta có công thức tường

minh cho xấp xỉ của cung trắc địa cho bởi (2.27) và giải một phương trình vi phân ma

trận khác (2.28). Xong lưu ý rằng (2.28) có 2k ẩn hàm (được xét như là phần tử trong

ma trận M (t) và N (t)) và là phương trình trình tuyến tính so với phương trình vi phân

(2.26), có nk hàm và là một phương trình phi tuyến.

Để ý rằng xấp xỉ cung trắc địa (để cho gọn, ta bỏ từ “xấp xỉ” trong những lần

Y0 ∈ V(k, n) và H ∈ TY0V(k, n). Ta sẽ ký hiệu nó bởi G(t, Y0, H). Trong trường hợp

k = 1, mêtric Euclide và mêtric chính tắc (2.25) trùng nhau. Khi đó đường trắc địa có

nhắc sau) nên trong Định lý 2.4.1 là xác định duy nhất khi biết điều kiện ban đầu

, t ≥ 0.

G(t, Y0, H) = cos (||H||t) Y0 + sin (||H||t)

dạng tường minh (xem thêm [2] để biết cách xây dựng).

H ||H||

(2.31)

2.4.3 Di chuyển song song

Trong không gian, việc tịnh tiến một véctơ rất dễ dàng vì véctơ thu được luôn là

một véctơ trong không gian tiếp xúc vốn trùng với không gian đang xét. Trên một đa

tạp không tầm thường, nó và không gian tiếp xúc của nó vốn là hai đối tượng khác

biệt. Việc di chuyển một véctơ tiếp xúc sao cho sau đó, nó vẫn là một véctơ tiếp xúc

30

đòi hỏi những thao tác đặc biệt. Thứ nhất, ta không thể tùy ý di chuyển véctơ mà phải

di chuyển dọc theo một cung trắc địa G(t) (vốn thay cho véctơ tịnh tiến trong không

gian). Tiếp đến ta phải loại bỏ các thành phần nằm trong không gian chuẩn tắc (tức

nằm ngoài không gian tiếp xúc). Cụ thể, ta cần giải phương trình di chuyển song song

˙∆ = −

trong đa tạp Stiefel

Γ (cid:0)∆, G, ˙G(cid:1) , ∆(0) = K ∈ TY0,

1 2

(2.32)

Γ(∆, G, ˙G) = (cid:0)∆ ˙GT + ˙G∆T (cid:1) G

với

+ G (cid:0) ˙GT (cid:0)In − GGT (cid:1) ∆ + ∆T (cid:0)In − GGT (cid:1) ˙G(cid:1)

(2.33)

Trong trường hợp tổng quát, ta không có công thức tường minh để biểu diễn nghiệm

của (2.32). Tuy nhiên, khi K = 1, di chuyển song song của K ∈ TY0 dọc theo cung

trắc địa G(t, Y0, H) được cho bởi

K.

∆(t) =

In + (cos(||H||t − 1))

HH (cid:62) ||H||2 − sin(||H||t)

Y0H (cid:62) ||H||

(cid:19) (cid:18) (2.34)

Ta sẽ ký hiệu véctơ đó là ∆(K, t, Y0, H).

2.4.4 Tích phân số

˙Y (t) = F (Y (t))Y (t), t > 0, Y0 ∈ V(k, n),

Xét phương trình vi phân trên V(k, n) có dạng

(2.35)

trong đó F : V(k, n) → Rn×n là một hàm ma trận bị chặn, Lipschitz và thỏa mãn điều

Y (cid:62) (cid:0)F (cid:62)(Y ) + F (y)(cid:1) Y = 0, ∀ Y ∈ V(k, n).

kiện phản đối xứng yếu

(2.36)

Từ Định lý 2.3.1 và cấu trúc (2.24) của không gian véctơ, phương trình (2.35) với điều

kiện (2.36) đảm bảo nó là một phương trình vi phân trên V(k, n).

31

Cũng như các phương pháp trình bày ở Chương 1, ta cũng cần chia khoảng thời

Yi ≈ Y (ti). Điều khác biệt ở đây là ta cần đảm bảo nếu Yi ∈ V(k, n) thì Yi+1 ∈ V(k, n).

gian thành những khoảng nhỏ (ti, ti+1) với ti+1 = ti + h, t0 = 0 và tính nghiệm xấp xỉ

F (Yi)Yi = ˙Yi,

Hiện tại, ta xét các phương pháp 1 bước: tính Yi+1 chỉ dựa trên Yi. Do điều kiện (2.36),

Yi + hY (cid:48)

i ta cập nhật theo cung trắc địa xuất phát từ Yi, hướng F (Yi)Yi và độ dài h.

cho ta hướng cập nhật. Trên đa tạp, thay cho phép cập nhật theo đường thẳng Yi+1 =

Ki+1 = F (Yi)Yi

Theo đó, mỗi bước giải phương trình vi phân trên V(k, n) gồm 2 bước phụ:

Yi+1 = G(h, Yi, Ki+1).

(2.37)

Khi k = 1, đường trắc địa không đòi hỏi tính toán vì đã có công thức tường minh

(2.31). Trong trường hợp tổng quát k > 1, ta cần giải phương trình (2.28) để có (2.27).

Thông thường để tăng tính chính xác, ta cần giải (2.28) bằng một phương pháp bậc

hai. Việc này sẽ không gây thêm khó khăn gì vì (2.28) chỉ là một phương trình vi phân

không ràng buộc, có thể được giải bằng bất kỳ một phương pháp bậc hai nào đã trình

bày trong Chương 1.

Tiếp theo, ta sẽ thảo luận phương pháp tích phân bậc hai. Những phương pháp này

chắc chắn sẽ cho kết quả tốt hơn phương pháp một bước vừa trình bày ở trên do nó kết

hợp được nhiều hướng tìm ý tưởng cụ thể như sau. Giả sử ta đang ở điểm Yi và hướng

tìm V(Yi), sử dụng phương pháp một bước trình bày ở trên, ta tìm được ˜Yi +1 và hướng

tìm V( ˜Yi + 1). Lẽ ra ta có thể sử dụng thông tin này để tìm nghiệm tại nút i + 2 theo

phương pháp một bước. Nhưng thay vì thế, ta di chuyển hướng tìm V( ˜Yi + 1), vốn là

một vectơ tiếp xúc tại ˜Yi + 1 về thành một vectơ tiếp xúc tại Yi bằng cách giải phương

trình dịch chuyển song song (2.32). Sau đó, kết hợp với V(Yi) đã biết, ta cập nhật lại

để tìm (Yi+1). Ta sẽ xem xét nó dưới góc độ của phương pháp kiểu Runge-Kutta bậc

hai.

32

2.4.5 Phương pháp kiểu Runge-Kutta bậc hai

0

c1

0

c2

a21

b1

b2

Xét phương pháp Runge-kutta kiểu 2 bước xác định bởi bảng Butcher sau đây

Sử dụng ý tưởng trình bày ở trên, phương pháp này áp dụng cho phương trình vi phân

trên đa tạp Stiefel có dạng

Algorithm 1 Phương pháp Runge-Kutta hai bước cho PTVP trên đa tạp

Stiefel 1: K1 = F (Y0)Y0;

2: Y01 = G (a21h, Y0, K1); 3: ˜K2 = F (Y01)Y01;

˙Y01 = ˙G (a21h, Y0, K1);

(cid:17)

hΓ

(cid:16) ˜K2, Y01, ˙Y01

;

4: 1 5: S2 = ˜K2 + a21 2 6: K2 = πτY0 (S2); 7: K = b1K1 + b2K2;

8: Y1 = G (h, Y0, K).

Ta có thể thấy trong Thuật toán 1, bước 2 và bước 4 lần lượt tính cung trắc địa sử

dụng công thức (2.27) và công thức đạo hàm của đường trắc địa (2.29). Các bước 5 (cid:1) thực chất là di chuyển song song của véctơ và 6 đỉnh xấp xỉ ∆ (cid:0)−a21h, ˜K2, Y01, ˙Y01 tiếp xúc tại Y0, ˜K2 dọc theo cung trắc địa G(t, Y01, ˜Y01). Phép chiếu 6 được đưa vào để

đảm bảo nghiệm của phương trình dịch chuyển song song (2.32) nằm trên không gian

tiếp xúc TY0Vn,k. Định lý sau đây chỉ ra độ chính xác bậc hai của Thuật toán 1.

K1 = F (Y0)Y0, Y1 = G (h, Y0, K1) ,

Định lí 2.4.2. Phương pháp kiểu Runge-Kutta

a21b2 = 1/2.

áp dụng vào hệ (2.35), với 1 < p ≤ n, có độ chính xác bậc hai nếu b1 + b2 = 1 và

33

Chứng minh. Bằng khai triển Taylor, nghiệm ∆ (cid:0)t, ˜K2, Y01, ˙Y01 (cid:1) của phương trình (2.32)

∆(0) = ˜K2 + a21

Γ (cid:0) ˜K2, Y01, ˙Y01

tại t = 0, có thể viết bởi

h 2

(cid:1) + O(h2) = S2 + O(h2).

K = b1K1 + b2K2 = b1K1 + b2∆(0) + O (cid:0)h2(cid:1)

Từ ∆(0) ∈ TY0 nên ∆(0) = πY0 [∆(0)] và ∆(0) = K2 + O(h2), vì vậy

= b1K1 + b2F (Y01) Y01 + b2a21

Γ (cid:0) ˜K2, Y01, ˙Y01

h 2

(2.38) (cid:1) + O (cid:0)h2(cid:1) .

Cho S(Y ) = F (Y )Y , do Y01 = Y0 + a21hK1 + O(h2), nên chúng ta có thể viết

S(Y01) = S(Y0) + a21hSY (K1) +

a2 21SY Y (K1) + · · · ,

h2 2

(2.39)

n (cid:88)

p (cid:88)

n (cid:88)

p (cid:88)

, SY Y =

SY (K1) =

(K1)ij(K1)ke.

∂S(Y0) ∂yij

∂2S(Y0) ∂yij∂yke

j,e=1

i=1

j=1

i=1,k

với

Ngoài ra, chúng ta giả sử

Γ (cid:0) ˜K2, Y01, ˙Y01

(cid:1) + O(h).

Do đó, bằng (2.33) chúng ta có

Γ (cid:0)K1, Y0, ˙Y0

1 Y0 + 2Y0K(cid:62) 1

0

(cid:1) = 2K1K(cid:62) (cid:0)In − Y0Y (cid:62) (cid:1) K1. (2.40)

Kết hợp (2.39) và (2.40), thì (2.38) dẫn đến

K = (b1 + b2) K1 + b2a21hSY (K1) (cid:2)K1K(cid:62)

+ b2ha21

0

1 Y0 + 2Y0K(cid:62) 1

(2.41) (cid:3) + O(h2). (cid:1) K1 (cid:0)In − Y0Y (cid:62)

Y1 = G (h, Y0, K) = Y0 + hK +

Bây giờ, từ (2.30) ta suy ra

0 K − Y0K(cid:62)K(cid:1) + O(h3),

0 K + KY (cid:62)

h2 2

(cid:0)Y0K(cid:62)Y0Y (cid:62)

và từ K = (b1 + b2) K1 + O(h) chúng ta có

Y1 = Y0 + hK +

(b1 + b2) (cid:2)Y0K(cid:62)

1 Y0Y (cid:62)

0 K1 + K1Y (cid:62)

0 K1 − y0K(cid:62)

1 K1

h2 2

(cid:3) + O(h3).

34

Y = Y0 + h (b1 + b2) K1 + b2a21h2SY (K1)

Vì vậy, do (2.41) nên

+b2a21h2 (cid:2)K1K(cid:62)

1 Y0 + Y0K(cid:62) 1

0

(cid:3) (cid:0)In − Y0Y (cid:62) (cid:1) K1

+

(b1 + b2) (cid:2)Y0K(cid:62)

1 K1

0 K1 − Y0K(cid:62)

0 K1 + K1Y (cid:62)

1 Y0Y (cid:62)

h2 2

= Y0 + h (b1 + b2) K1 + b − 2a21h2SY (K1)

+

[(b1 + b2 − 2b2a21)YoK(cid:62)

1 K1

0 K1 + (2b2a21 − b1 − b2)Y0K(cid:62)

1 Y0Y (cid:62)

h2 2 +K1Y (cid:62)

0 (2b2a21F (cid:62)(Y0) + (b1 + b2)F (Y0))Y0] + O(h3).

(cid:3) + O(h3)

Y (h) = Y0 + h ˙Y0 +

SY ( ˙Y0) + O(h3).

h2 2

Vì vậy, do F (Y ) thỏa mãn điều kiện đối xứng lệch yếu (2.36) chúng ta có

Vì vậy nếu b1 + b2 = 1 và 2b2a21 = 1 ta có điều cần chứng minh.

Chú ý 2.4.3. Ngay cả khi ta áp dụng một lược đồ bậc cao ta cũng không thể tăng độ

chính xác kể trên. Lí do là vì khi giải xấp xỉ phương trình dịch chuyển song song ta

chỉ sử dụng phương pháp có độ chính xác bậc hai. Trong mục sau, ta sẽ xét trường

hợp k = 1 khi đa tạp Sitefel trở thành mặt cầu đơn vị.

2.4.6 Phương pháp giải số trên mặt cầu đơn vị

Điều khác biệt so với trường hợp tổng quát là ta có công thức hiện cho cung trắc

địa và dịch chuyển song song. Xét phương pháp Runge-Kutta hiện s bước định nghĩa

c A

b

bởi bảng Butcher

trong đó c = (ci), b = (bi) và A = (aij) với các phần tử 0 trừ các phần tử hàng i + 1,

cột i, i = 1, · · · , s − 1. Khi đó ta có thuật toán.

35

Algorithm 2 Phương pháp Runge-Kutta giải số trên mặt cầu đơn vị 1: K1 = F (Y0) Y0; 2: for i = 1, · · · , s − 1 do

3:

4:

5:

(cid:17)

6:

Y0i = G (ai+1,ih, Y0, Ki); ˜Ki+1 = F (Y0i) Y0i; ˙Y0i = ˙G (ai+1,ih, Y0, Ki); Ki+1 = ∆

(cid:16) ˜Ki+1, −ai+1,ih, Y0i, ˙Y0i

;

7: end for s (cid:88)

8: K =

biKi;

i=1 9: Y1 = G (h, Y0, K).

Ta lưu ý rằng cung trắc địa ở 3 và dịch chuyển song song ở bước 6 sử dụng các công

thức hiện (2.31) và (2.34). Định lý sau đây cho ta thông tin về tính chính xác của thuật

toán.

Định lí 2.4.4. Thuật toán 2 áp dụng cho (2.35) với p = 1 bảo toàn bậc chính xác của

sơ đồ Runge-Kutta cơ bản.

Chứng minh. Sử dụng kí hiệu của Định lí (2.4.2), khai triển Taylor của điều kiện bậc

Y1 = Y0 + h (b1 + b2) K1 + h2b2a21SY (K1) + h2 (b2a21 − 1/2) (cid:107)K1(cid:107)2 Y0 + O (cid:0)h3(cid:1) .

hai của lược đồ Runge-Kutta bậc hai- giai đoạn (5) là

Bằng điều kiện bậc hai của lược đồ Runge-Kutta bậc hai nên Y1 = Y (h) = O(h3).

Sử dụng tính chất phép đối xứng của tích vô hướng trên Sn, khai triển Taylor của

36

3 (cid:88)

Y1 =Y0 + (

bi)K1 + h2 (b2a21 + b3a32) (K1)

i=1 

nghiệm của phương pháp Runge-Kutta bậc 3 giai đoạn 5 là

+ h2

bi

1 2

i=1

(cid:33)2 (cid:32) 3 (cid:88)  ||K1||2Y0 b2a21 + b3a32 −

+

bi

32 −

21 + b3a2 32

21 + b3a2

1 2

h3 2

h3 2

i=1 (cid:35)

(cid:33)3 (cid:32) 3 (cid:88) (cid:0)b2a2 (cid:1) SY Y (K1) +  ||K1||2K1  b2a2

+ h3

(b2a21 + b3a32)

b2a2

bi

SY (K1)(cid:62)K1Y0

21 + b3a2

32 + b3a32a21 −

i=1

+ O(h4).

(cid:33) (cid:34) (cid:32) 3 (cid:88)

Kỹ thuật tương tự có thể được dùng để chứng minh bậc chính xác p bất kỳ của lược

đồ s bước.

2.5 Phương trình vi phân trên nhóm Lie

Trong mục này, chúng tôi trình bày một số vấn đề liên quan đến nhóm Lie như: sơ

lược về nhóm Lie, phương trình vi phân trên nhóm ma trận Lie,. . .. Những nội dung

này chủ yếu dựa vào tài liệu [5].

2.5.1 Sơ lược về nhóm Lie

Một nhóm Lie G là một đa tạp khả vi và đồng thời là một nhóm mà phép toán

trên đó, tạm gọi là phép nhân, là một ánh xạ khả vi.

Nhóm Lie là khái niệm tương đối rộng, song trong khuôn khổ của luận văn này, ta

chỉ xét các nhóm Lie ma trận, tức là các nhóm Lie mà là nhóm con của GL(n), nhóm

các ma trận khả nghịch với phép nhân ma trận thông thường.

O(n) = {Y ∈ GL(n) : Y (cid:62)Y = Y Y (cid:62) = I}

Ví dụ 2.5.1. Nhóm Lie các ma trận trực chuẩn

37

Ta có thể dễ dàng chứng minh được đây là một đa tạp con nhúng trong Rn×n bằng

việc xét ánh xạ Y (cid:55)−→ Y (cid:62)Y − I từ Rn×n vào sym(n) và áp dụng Định lý hàm ẩn trong

hình học vi phân. Hơn thế nữa phép nhân ma trận là phép toán cơ bản nên khả vi. Do

Bảng 2.1: Một vài nhóm ma trận Lie & đại số Lie tương ứng

đó O(n) là một nhóm Lie với phần tử đơn vị là ma trận đơn vị I.

Nhóm Lie Đại số Lie

GL(n)= {Y | det Y (cid:54)= 0} gl(n)= {A| ma trận bất kỳ}

nhóm đơn tuyến tính tổng quát nhóm đại số của ma trận m × n

SL(n)= {Y | det Y = 1} sl(n)= {A| vết (A)}

nhóm tuyến tính đặc biệt đại số Lie tuyến tính đặc biệt

so(n)= (cid:8)A| AT + A = 0(cid:9) OL(n)= (cid:8)Y | Y T Y = I(cid:9)

nhóm trực giao ma trận đối xứng lệch

so(n)= (cid:8)A| AT + A = 0(cid:9) SO(n)= {Y ∈ O(n)| detY = 1}

nhóm trực giao đặc biệt

ma trận đối xứng lệch sp(n)= (cid:8)A| JA + AT J = 0(cid:9) Sp(n) = (cid:8)Y | Y T JY = J(cid:9)

nhóm ngẫu đối

Ký hiệu g = TI G là không gian tiếp xúc của một nhóm ma trận Lie G tại phần tử đơn

vị. Ta có mệnh đề sau.

[A, B] = AB − BA,

Mệnh đề 2.5.2. Móc Lie trên tập g = TI G, định nghĩa là

xác định một toán tử trong g, song tuyến tính, phản đối xứng và thỏa mãn đẳng thức

[A, [B, C]] + [C, [A, B]] + [B, [C, A]] = 0.

Jacobi

Do những kết luận của Mệnh đề 2.5.2, g trở thành một đại số và được gọi là đại số Lie

của nhóm Lie G. Trước khi kết thúc mục này ta phát biểu mệnh đề sau.

38

exp : g → G

A (cid:55)→ eA.

Mệnh đề 2.5.3. Ánh xạ mũ ma trận

là một vi phôi địa phương trong một lân cận của A = 0.

2.5.2 Phương trình vi phân trên nhóm ma trận Lie

Cho G là một nhóm ma trận Lie, giả sử Y ∈ G bất kì, A ∈ g. Khi đó, theo định

˙γ(0) = A. Xét cung

γ(t) := α(t)Y.

nghĩa, phải tồn tại cung trơn α(t) trên G sao cho α(0) = I và

˙γ(0) = ˙α(0)Y = AY . Như vậy, không gian tiếp xúc của nhóm

G tại một điểm Y bất kỳ có dạng

TY G = {AY : A ∈ g}.

Rõ ràng γ(t) ∈ G, ∀t và

˙Y = A(Y )Y.

Bây giờ, ta hãy xét phương trình vi phân dạng

(2.42)

Theo Định lý 2.3.1, phương trình (2.42) là một phương trình vi phân trên G nếu nó

thỏa mãn A(Y )Y ∈ TY G hay A(Y ) ∈ g.

G = {Y : g(Y ) = const},

Nếu ngoài điều kiện A(Y ) ∈ g với mọi y ta còn thêm

là một nhóm Lie ma trận thì g(Y ) sẽ trở thành một bất biến của phương trình (2.42)

theo Định nghĩa ở 2.5.1.

2.5.3 Phương pháp Crouch-Grossmann

˙Y = A(Y )Y, Y (0) = Y0,

Trong mục này, ta sẽ xét phương trình vi phân

(2.43)

39

thỏa mãn Y0 thuộc nhóm Lie ma trận G và A(Y ) là một phần tử của đại số Lie g tương

ứng của G. Đương nhiên ta có thể sử dụng phương pháp chiếu như trình bày ở Mục

2.3 để giải (2.43). Tuy nhiên, ở đây ta sẽ xét một phương pháp nội tại, tức cho ta xấp

xỉ trên đa tạp đang xét.

Như đã biết, phương pháp Runge-Kutta hiện về cơ bản gồm hai thao tác chính

(i) Tính trường vectơ f (Y ) = A(Y )Y .

(ii) Cập nhật nghiệm ở dạng Y + haf (Z).

Ta có thể quan sát thấy trong trường hợp phương trình vi phân trên nhóm Lie, ngay cả

khi Y và Z ∈ G, cập nhật Y + haA(Z)Z có thể không thuộc G. Ý tưởng của Crouch

và Grossmann là thay thế cập nhật đó bởi exp (haA(X)) Y .

Định nghĩa 2.5.4. Cho bi, aij, i, j = 1, · · · s là các số thực. Một phương pháp Crough-

Grossmann s-bước được cho bởi

Y (i)(cid:17)

,

Y (i) = exp (hai,i−1Ki−1) · . . . · exp (hai1K1) Yn, Ki = A

(cid:16)

Yn+1 = exp (hbsks) · . . . · exp (hb1k1) Yn.

(2.44)

Chẳng hạn với s = 2, a21 = 1, b1 = b2 = 1/2, phương pháp trong Định nghĩa 2.5.4 là

exp

K2

K1

Yn,

Yn+1 = exp

(cid:19) (cid:19)

(cid:18)h 2 (cid:18)h 2

với K1 = A(Yn), K2 = A (exp(hki)Yn).

Rõ ràng do các nhân tử trong (2.44) đều thuộc G nên phương pháp này cho một

xấp xỉ nghiệm luôn nằm trên G. Ta hãy xem phương pháp này chính xác như thế nào.

j aij. Một phương pháp Crouch-Grossman có bậc p

Định lí 2.5.5. Giả sử ci = (cid:80)

(p ≤ 3) nếu các điều kiện tương ứng sau được thỏa mãn:

bi = 1

Bậc 1: (2.45)

bici = 1/2

Bậc 2: (2.46) (cid:88) i (cid:88) i

40

bic2

3 = 2/3

Bậc 3: (2.47)

biaijcj = 1/6

ij

(cid:88) i (cid:88) (2.48)

bicibj = 1/3.

b2 i ci + 2

i

(cid:88) (2.49) (cid:88) i

Chứng minh. Như trường hợp phương pháp Runge-Kutta, điều kiện bậc có thể được

tìm thấy bằng cách so sánh khai triển chuỗi Taylor chính xác và lời giải bằng số. Ngoài

điều kiện được phát biểu trong định lý này, nó còn dẫn đến quan hệ như là

.

bibjcj =

b2 i ci + 2

2 3

i

i

(cid:88) (cid:88) (2.50)

ij bicibj = 1, thỏa

Thêm phương trình này vào phương trình (2.49) chúng ta tìm 2 (cid:80)

mãn bởi (2.45) và (2.46). Vì thế, quan hệ (2.50) là hệ quả của điều kiện được phát

biểu trong định lý.

Để tiện sử dụng, ta có cụ thể hóa các điều kiện trong Định lý 2.5.5 thành bảng

Bảng 2.2: Bảng 1

Bảng 2.3: Bảng 2

0

0

−1/24 −1/24

3/24

3/4

161/24 −6

17/24

17/24 119/216 17/108

1

−2/3 2/3

13/51

−2/3

24/17

sau:

2.6 Ví dụ số

Trong mục này, chúng tôi trình bày một số ví dụ để minh họa cho một vài khía

cạnh lý thuyết của luận văn.

41

Trước tiên, ta xem xét mô hình của vật rắn không biến dạng đã được nhắc đến

˙y1 = a1y2y3, a1 = (I2 − I3) / (I2I3)

trong Ví dụ 2.1.7.

˙y2 = a2y3y1, a2 = (I3 − I1) / (I3I1)

˙y3 = a3y1y2, a3 = (I1 − I2) / (I1I2) ,

(2.51)

trong đó, y = (y1, y2, y3)(cid:62) là mômen động lượng (argular momentum) và I1, I2, I3 là

các mômen quán tính.

H1 (y1, y2, y3) = y2

Như đã phân tích, hệ cơ học này có bất biến

1 + y2

2 + y2

3 − C1,

(2.52)

−

+

+

C2.

H2 (y1, y2, y3) =

1 2

1 2

y2 2 I2

y2 3 I3

√

C1

và (cid:19) (2.53) (cid:18) y2 1 I1

√

√

I1C2,

I2C2 và

Khi cho C1, C2 các giá trị dương cụ thể, ta thấy (2.52) mô tả mặt cầu bán kính

I3C2. Như vậy, về mặt hình học, các đường cong tích phân của (2.51) phải đồng thời

trong khi (2.53) là một ellipse có tâm tại gốc tọa độ và có bán trục √

nằm trên một hình cầu và một ellipse. Nói cách khác, nó chính là giao tuyến của hai

mặt đó.

Để giải số phương trình vi phân (2.51), chúng tôi chọn thời điểm ban đầu t0 = 1

với bước thời gian h = 0.13.

42

Hình 2.1: Sử dụng phương pháp Euler tiến và không chiếu.

Trong Hình 2.1, chúng tôi minh họa nghiệm thu được khi sử dụng phương pháp Euler

tiến. Có thể dễ dàng quan sát thấy sau một vài bước, nghiệm số này đã nhảy ra khỏi

mặt cầu. Điều này một lần nữa cảnh báo việc giải phương trình vi phân mà không

xem xét các bất biến của nó cũng như không sử dụng các phương pháp hữu hiệu để

Hình 2.2: Sử dụng phương pháp Euler tiến và kết hợp.

bảo toàn các bất biến đó.

Trong Hình 2.2, chúng tôi giữ nguyên các thông số như hình trên nhưng sử dụng

43

phương pháp chiếu trình bày ở Mục 2.3 để tìm nghiệm số và minh họa nghiệm số đó

bằng màu xanh. Do đa tạp là mặt cầu nên phép chiếu đơn giản chỉ là giao của đường

thẳng nối tâm và điểm biểu diễn nghiệm với mặt cầu.

[1, 1, 1](cid:62) và A là

V(1, 3). Để xây dựng nghiệm giải tích, chúng tôi đặt t0 = 0, q0 =

1 √ 3

Ở ví dụ thứ hai, chúng tôi xét mặt cầu đơn vị như là một mô hình cho đa tạp Stiefel

ma trận đường chéo có các phần tử chéo là a1 = −α ; a2 = α; a3 = −α với α là một

˙q(t) = Aq(t),

hằng số dương. Khi đó, dễ thấy phương trình vi phân

(2.54)

e−αt

q(t) =

.

có nghiệm  

eαt

1 √ 3

(2.55)

e−αt

         

Để thu được các nghiệm trên mặt cầu đơn vị, ta chuẩn hóa Y (t) = q(t)/||q(t)||. Khi

˙Y = F (Y )Y = (cid:0)I3 − Y Y (cid:62)(cid:1) AY, t > 0, Y (0) = q0.

đó, có thể thấy Y (t) chính là nghiệm của phương trình vi phân

(2.56)

Vậy, ta đã thu được một phương trình vi phân trên mặt cầu đơn vị (2.56).

Ta giải số phương trình này trên khoảng [0; 10] với bước thời gian 0.01. Ta sử dụng

phương pháp đã trình bày ở Mục 2.4.4 có bậc 1 để giải. Nghiệm chính xác của nó (tô

màu đỏ) được so sánh với nghiệm chính xác (tô màu xanh) trong Hình 2.3.

44

Hình 2.3: So sánh nghiệm số và nghiệm chính xác trên đa tạp Stiefel V(1, 3).

Ta có thể quan sát thấy bên cạnh việc bảo toàn bất biến (nằm trên mặt cầu đơn

vị), nó gần như trùng với nghiệm chính xác. Người đọc có thể tham khảo sai số tuyệt

Hình 2.4: Sai số tuyệt đối của nghiệm số trên đa tạp Stiefel V(1, 3).

đối so với nghiệm chính xác trình bày trong Hình 2.4.

45

Kết luận

Trong luận văn này, chúng tôi đã trình một số phương pháp số thông dụng giải

phương trình vi phân ma trận với ràng buộc đa tạp. Cụ thể, chúng tôi trình bày khái

niệm bất biến tổng quát, bất biến bậc một và bất biến bậc hai. Sau đó, chúng tôi làm

rõ với điều kiện gì thì các phương pháp tích phân phương trình vi phân thông thường

như phương pháp Runge-Kutta, phương pháp Nystr¨om bảo toàn các bất biến này. Tiếp

đó, chúng tôi trình bày những kỹ thuật cụ thể để tích phân phương trình vi phân trên

hai đạp Stiefel và nhóm Lie. Cuối cùng, chúng tôi xét một số ví dụ số để minh họa

cho phương pháp được trình bày.

46

Tài liệu tham khảo

Tiếng Việt

[1] Khu Quốc Anh, Nguyễn Doãn Tuấn (2005), Lí thuyết liên thông và hình học

Riemann, NXB Đại học Sư Phạm.

[2] Hoàng Ngọc Thế (2009), Về tối ưu trên đa tạp Riemann. Luận văn thạc sĩ toán

học, Đại học Khoa học, Đại học Thái Nguyên.

Tiếng Anh

[3] L.D. Buono and L. Lopez (2001), "Runge-Kutta type methods based on

geodesics for systems of odes on the Stiefel manifold", BIT numerical mathe-

matics. Appl, 41(5): 912–923.

[4] A. Edelman, T.A. Arias, and S.T. Smith (1998), "The geometry of algorithms

with orthogonality constraints", SIAM J. Matrix Anal. Appl., 20(2), 303–353.

[5] E. Hairer, G. Wanner and C. Lubich (2001), Geometric numerical integration:

structure-preserving algorithms for ordinary differential equations, Springer.