Đăng ký Đăng nhập
Trang chủ Tìm nghiệm phương trình đạo hàm riêng bằng phương trình sai phân...

Tài liệu Tìm nghiệm phương trình đạo hàm riêng bằng phương trình sai phân

.PDF
63
41
64

Mô tả:

ĐẠI HỌC THÁI NGUYÊN TRƯỜNG ĐẠI HỌC KHOA HỌC ------------------------------- NGUYỄN DUY TIN TÌM NGHIỆM PHƯƠNG TRÌNH ĐẠO HÀM RIÊNG BẰNG PHƯƠNG TRÌNH SAI PHÂN 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 ------------------------------- NGUYỄN DUY TIN TÌM NGHIỆM PHƯƠNG TRÌNH ĐẠO HÀM RIÊNG BẰNG PHƯƠNG TRÌNH SAI PHÂN 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 THỊ NGỌC OANH THÁI NGUYÊN - 2019 1 Mục lục Trang Danh sách hình vẽ 2 Lời nói đầu 3 Chương 1 Một số kiến thức cơ bản 5 1.1. Toán tử sai phân . . . . . . . . . . . . . . . . . . . . . . . 5 1.2. Tính tổng . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3. Biến đổi z . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Chương 2 Giải phương trình đạo hàm riêng bằng phương trình sai phân 19 2.1. Rời rạc hóa phương trình đạo hàm riêng . . . . . . . . . . 19 2.2. Nghiệm của phương trình đạo hàm riêng . . . . . . . . . . 27 2.3. Ví dụ số . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.1. Phương trình parabolic tuyến tính 1 chiều . . . . . 35 2.3.2. Phương trình parabolic tuyến tính 2 chiều . . . . . 39 Phụ lục 60 2 Danh sách hình vẽ 2.1 Mô hình đối với phương trình truyền nhiệt . . . . . . . . . 21 2.2 Lưới điểm đạt được từ trục ban đầu . . . . . . . . . . . . . 22 2.3 Lưới điểm đạt được không với giá trị biên . . . . . . . . . 22 2.4 Mô hình phương pháp ẩn . . . . . . . . . . . . . . . . . . . 25 2.5 Mô hình phương trình Laplace . . . . . . . . . . . . . . . . 27 2.6 Nghiệm xấp xỉ Ví dụ 2.5 . . . . . . . . . . . . . . . . . . . 37 2.7 Nghiệm chính xác Ví dụ 2.5 . . . . . . . . . . . . . . . . . 38 2.8 Sai số Ví dụ 2.5 . . . . . . . . . . . . . . . . . . . . . . . . 38 2.9 Nghiệm xấp xỉ Ví dụ 2.6 . . . . . . . . . . . . . . . . . . . 39 2.10 Nghiệm chính xác Ví dụ 2.6 . . . . . . . . . . . . . . . . . 40 2.11 Sai số Ví dụ 2.6 . . . . . . . . . . . . . . . . . . . . . . . . 40 2.12 Nghiệm xấp xỉ Ví dụ 2.7 . . . . . . . . . . . . . . . . . . . 43 2.13 Nghiệm chính xác Ví dụ 2.7 . . . . . . . . . . . . . . . . . 43 2.14 Sai số Ví dụ 2.7 . . . . . . . . . . . . . . . . . . . . . . . 44 3 Lời nói đầu Sử dụng phương trình sai phân để giải số là một phương pháp khá phổ biến và hữu hiệu khi nghiên cứu mô hình toán học liên quan tới các vấn đề trong khoa học, kỹ thuật, kinh tế và nhiễu lĩnh vực khác của thực tiễn, ta có thể tìm thấy rất nhiều ví dụ cụ thể như trong Chương 1 tài liệu [2]. Phương pháp này được đề xuất từ nửa cuối những năm 40 của thế kỷ trước và ngày càng khẳng định vai trò quan trọng trong giải tích ứng dụng và đặc biệt quan trọng trong việc nghiên cứu nghiệm số của phương trình đạo hàm riêng [1, 3]. Luận văn trình bày một số kiến thức cơ bản liên quan tới phương trình sai phân và áp dụng phương trình sai phân tìm nghiệm số phương trình đạo hàm riêng tuyến tính. Luận văn được chia làm hai chương. Chương 1 là chương mở đầu, trình bày các kiến thức cơ bản nhất liên quan tới phương trình sai phân như: Định nghĩa và các tính chất của toán tử sai phân; định nghĩa và các tính chất của tổng bất định; biến đổi z giải phương trình sai phân. Nội dung của Chương 1 được tham khảo chủ yếu trong hai tài liệu [1, 2]. Chương 2 nghiên cứu ứng dụng phương trình sai phân vào việc tìm nghiệm số của phương trình đạo hàm riêng. Chương này cũng trình bày các bước rời rạc bài toán, nghiệm rời rạc đồng thời có các ví dụ số minh họa thông qua ngôn ngữ lập trình MATLAB. Nội dung của chương này được tham khảo chủ yếu trong hai tài liệu [1, 3]. Sau thời gian học tập và rèn luyện tại Trường Đại học Khoa học – 4 Đại học Thái Nguyên, bằng sự biết ơn và kính trọng, em xin gửi lời cảm ơn chân thành đến Ban Giám hiệu, Phòng Đào tạo, Khoa Toán - Tin Trường Đại học Khoa học – Đại học Thái Nguyên và các thầy cô đã nhiệt tình hướng dẫn, giảng dạy và tạo mọi điều kiện thuận lợi giúp đỡ em trong suốt quá trình học tập, nghiên cứu và hoàn thiện đề tài luận văn Thạc sĩ này. Đặc biệt, em xin bày tỏ lòng biết ơn sâu sắc tới Tiến sĩ Nguyễn Thị Ngọc Oanh, người đã trực tiếp hướng dẫn, giúp đỡ em trong quá trình thực hiện đề tài. Xin chân thành cảm ơn gia đình, bạn bè cùng đồng nghiệp đã tạo điều kiện giúp đỡ, động viên trong quá nghiên cứu và hoàn thiện đề tài này. Em xin trân trọng cảm ơn! Thái Nguyên, ngày 01 tháng 11 năm 2019 Học viên Nguyễn Duy Tin 5 Chương 1 Một số kiến thức cơ bản Trong chương này chúng tôi trình bày một số kiến thức cơ bản liên quan tới định nghĩa và tính chất toán tử sai phân, tổng bất định và biến đổi z. 1.1. Toán tử sai phân Định nghĩa 1.1 Cho y(t) là một hàm biến thực t, toán tử sai phân ∆ được xác định bởi ∆y(t) = y(t + 1) − y(t). Ví dụ 1.1 Phương pháp Euler xấp xỉ nghiệm của bài toán ban đầu 0 x (t) = f (t, x(t)), (1.1) x(t0 ) = x0 , (1.2) 0 đạt được bằng cách thay x (t) bởi x(t + h) − x(t) . Ta có h x(t + h) − x(t) = f (t, x(t)) h hay x(t + h) = x(t) + hf (t, x(t)). Ta sẽ sử dụng một dạng thuận tiện hơn cho phương trình sai phân này, đặt xn = x(t0 + nh), n = 0, 1, 2, . . . , khi đó phương trình trên có thể 6 được viết lại dưới dạng sau xn+1 = xn + hf (t, xn (t)), n = 0, 1, 2, . . . trong đó x0 đã biết. Xét toán tử sai phân với bước lưới h > 0 là z(s + h) − z(s). Đặt y(t) = z(th), khi đó z(s + h) − z(s) = z(th + h) − z(th) = y(t + 1) − y(t) = ∆y(t). Sai phân cấp cao hơn được định nghĩa bằng hợp thành của toán tử sai phân và chính nó. Sai phân cấp hai được định nghĩa như sau ∆2 y(t) = ∆(∆y(t)) = ∆(y(t + 1) − y(t)) = (y(t + 2) − y(t + 1)) − (y(t + 1) − y(t)) = y(t + 2) − 2y(t + 1) + y(t). Sai phân cấp n được xây dựng theo công thức quy nạp như sau n(n − 1) ∆n y(t) = y(t + n) − ny(t + n − 1) + y(t + n − 2) + · · · + (−1)n y(t) 2! ! n X n = (−1)k y(t + n − k). (1.3) k k=0 Định nghĩa 1.2 Toán tử "dịch chuyển" được định nghĩa như sau Ey(t) = y(t + 1). Nếu I là toán tử đồng nhất, tức là Iy(t) = y(t) thì ta có ∆ = E − I. Các tính chất cơ bản của toán tử ∆ được cho bởi định lý sau đây (1.4) 7 Định lý 1.1 Toán tử sai phân ∆ thỏa mãn các tính chất sau đây a. ∆m (∆n y(t)) = ∆m+n y(t) với mọi số nguyên dương m và n. b. ∆(y(t) + z(t)) = ∆y(t) + ∆z(t). c. ∆(Cy(t)) = C∆y(t). d. ∆(y(t)z(t)) = y(t)∆z(t) + Ez(t)∆y(t).  y(t)  z(t)∆y(t) − y(t)∆z(t) = . e. ∆ z(t) z(t)Ez(t) Chứng minh. a. Ta chứng minh quy nạp theo m, đẳng thức hiển nhiên đúng với m = 1. Giả sử đẳng thức đúng với m = k, ta chứng minh đúng tới m = k + 1. Thật vậy ∆k+1 (∆n y(t)) = ∆k (∆(∆n y(t))) = ∆k (∆n+1 y(t))) = ∆k+n+1 y(t). b. Ta có ∆(y(t) + z(t)) = y(t + 1) + z(t + 1) − y(t) − z(t) = y(t + 1) − y(t) + z(t + 1) − z(t) = ∆y(t) + ∆z(t). c. Ta có ∆(Cy(t)) = Cy(t + 1) − Cy(t)  = C y(t + 1) − y(t) = C∆y(t). d. Ta có ∆(y(t)z(t)) = y(t + 1)z(t + 1) − y(t)z(t) = y(t + 1)z(t + 1) − y(t)z(t + 1) + y(t)z(t + 1) − y(t)z(t) = z(t + 1)∆y(t) + y(t)∆z(t) = y(t)∆z(t) + Ez(t)∆y(t). 8 e. Ta có  y(t)  y(t + 1) y(t) ∆ = − z(t) z(t + 1) z(t) y(t + 1)z(t) − z(t + 1)y(t) = z(t)z(t + 1)  y(t + 1)z(t) − y(t)z(t) − z(t + 1)y(t) − y(t)z(t) = z(t)Ez(t) z(t)∆y(t) − y(t)∆z(t) . = z(t)Ez(t) Định lý 1.2 Cho a là một hằng số. Khi đó a. ∆at = (a − 1)at . b. ∆ sin at = 2 sin a2 cos a(t + 12 ). c. ∆ cos at = −2 sin a2 sin a(t + 12 ). d. ∆ log at = log(1 + 1t ). e. ∆ log Γ(t) = log t, trong đó Γ(t) được xác định như sau Z ∞ e−z z t−1 dz. Γ(t) = 0 Chú ý 1.1 Hàm Γ(t) thỏa mãn tính chất Γ(t + 1) = t. Γ(t) Thật vậy, ta có Z ∞ e−z z t dz Γ(t + 1) = 0 = [−e−z z t ]∞ 0 Z =t ∞ Z − ∞ (−e−z )tz t−1 dz 0 e−z z t−1 dz = tΓ(t). 0 Chứng minh. a. Ta có ∆at = at+1 − at = (a − 1)at . 9 b. Ta có ∆ sin at = sin a(t + 1) − sin at a 1 = 2 sin cos a(t + ). 2 2 c. Ta có ∆ cos at = cos a(t + 1) − cos at 1 a = −2 sin sin a(t + ). 2 2 d. Ta có ∆ log at = log a(t + 1) − log at = log 1 a(t + 1) = log(1 + ). at t e. Ta có ∆ log Γ(t) = Γ(t + 1) − Γ(t) = log Γ(t + 1) Γ(t) = log t. Định nghĩa 1.3 Bậc giai thừa lùi tr theo giá trị r (đọc là t lùi về r) được xác định như sau: a. Nếu r = 1, 2, 3, . . . , thì tr = t(t − 1)(t − 2) · · · (t − r + 1). b. Nếu r = 0 thì t0 = 1. 1 . (t + 1)(t + 2) · · · (t − r) d. Nếu r không phải là số nguyên dương thì c. Nếu r = −1, −2, −3, . . . , thì tr = tr = Γ(t + 1) Γ(t − r + 1) Ta chú ý rằng, Định nghĩa 1.3 của tr chỉ xác định khi các giá trị t và r làm cho các công thức có nghĩa. Định nghĩa 1.4 Hệ số nhị thức t r ! t r ! được xác định bởi tr = . Γ(r + 1) 10 Hệ số nhị thức có tính chất sau đây • Tính đối xứng t ! r • Chuyển ra ngoài dấu ngoặc ! t r = = t r t ! . t−r t−1 r−1 ! . • Công thức cộng t r ! = t−1 r ! + t−1 ! r−1 . Định lý dưới đây sẽ chỉ ra một số kết quả liên quan tới sai phân. Định lý 1.3 a. ∆t tr =!rtr−1 . ! t t b. ∆t = , r 6= 0. r r−1 ! ! r+t r+t c. ∆t = . t t+1 Chứng minh. a. Trước tiên ta chứng minh cho trường hợp r là nguyên dương. ∆t tr = (t + 1)r − tr = (t + 1)t · · · (t − r + 2) − t(t − 1) · · · (t − r + 1) = t(t − 1) · · · (t − r + 2) [(t + 1) − (t − r + 1)] = rtr−1 . Trong trường hợp r là bất kỳ, sử dụng công thức d trong Định nghĩa 11 1.3, ta có ∆t tr = (t + 1)r − tr Γ(t + 2) Γ(t + 1) − Γ(t − r + 2) Γ(t − r + 1) (t + 1)Γ(t + 1) (t − r + 1)Γ(t + 1) = − Γ(t − r + 2) Γ(t − r + 2) Γ(t + 1) = rtr−1 . =r Γ(t − r + 2) = b. Ta có ∆t t ! r tr rtr−1 = = ∆t Γ(r + 1) Γ(r + 1) ! r−1 t t = = . Γ(r) r−1 c. Ta có ∆t r+t t ! = ∆t (r + t)t (r + t + 1)t+1 (r + t)t = − Γ(t + 1) Γ(t + 2) Γ(t + 1) Γ(r + t + 1) Γ(r + t + 2) − Γ(r + 1)Γ(t + 2) Γ(r + 1)Γ(t + 1) Γ(r + t + 2) − (t + 1)Γ(r + t + 1) = (t + 1)Γ(r + 1)Γ(t + 1) rΓ(r + t + 1) = Γ(r + 1)Γ(t + 2) ! t+1 r+t Γ(r + t + 1) (r + t) = = = . Γ(r)Γ(t + 2) Γ(t + 2) t+1 = Ví dụ 1.2 Tìm nghiệm của phương trình sai phân sau y(t + 2) − 2y(t + 1) + y(t) = t(t − 1). Phương trình sai phân trên có thể viết gọn dưới dạng ∆2 y(t) = t2 . Mặt khác, áp dụng công thức a của Định lý 1.3, ta có ∆2 t4 = ∆4t3 = 12t2 . 12 Do đó 4 ∆2 t4 2t t = =∆ = ∆2 y(t). 12 12 2 Do vậy t4 y(t) = 12 là một nghiệm của phương trình sai phân đã cho. 1.2. Tính tổng Định nghĩa 1.5 Tổng bất định (hay đối sai phân - antidifference) của P hàm y(t) được ký hiệu là y(t), là một hàm sao cho X  ∆ y(t) = y(t) thỏa mãn với mọi t trong miền xác định của hàm y. Định lý 1.4 Nếu z(t) là một tổng bất định của y(t), thì mọi tổng bất định của y(t) được cho bởi X y(t) = z(t) + C(t) trong đó C(t) có miền xác định giống z(t) và ∆C(t) = 0. Chứng minh. Rõ ràng, nếu z(t) là một tổng bất định của y(t), và ∆(t) = 0 ta có ∆(z(t) + C(t)) = ∆z(t) + ∆C(t) = y(t). Ngược lại, nếu f (t) là một tổng bất định của y(t) thì   ∆ f (t) − z(t) = ∆f (t) − ∆z(t) = 0. Do đó f (t) = z(t) + C(t), với ∆C(t) = 0. Ví dụ 1.3 Tính tổng S = P 6t . Từ Tính chất a của Định lý 1.2 ta có 6t ∆6 = 5.6 và ∆ = 6t . 5 t t 13 6t Do đó là một tổng bất định của 6t . Như vậy tất cả tổng bất định của 5 6t t 6 có dạng + C(t). Ta sẽ đi xác định hàm C(t). 5 Xét trường hợp đầu tiên khi miền xác định của y(t) là tập các số tự nhiên, khi đó ta có ∆C(t) = C(t + 1) − C(t) = 0. Từ đó ta có C(1) = C(2) = · · · , tức là C(t) = C với C là hằng số. Xét trường hợp miền xác định của y(t) là tập tất cả các số thực, khi đó từ phương trình ∆C(t) = C(t + 1) − C(t) = 0 ta rút ra C(t + 1) = C(t) với mọi số thực t hay C(t) là hàm tuần hoàn chu kì 1. Hệ quả 1.1 Cho y(t) là hàm xác định trên tập có dạng {a, a + 1, a + 2, . . . } trong đó a là số thực bất kì. Giả sử z(t) là một tổng bất định của y(t). Khi đó mọi tổng bất định của y(t) có dạng X y(t) = z(t) + C với C là hằng số bất kỳ. Định lý 1.5 Cho a là một hằng số và ∆C(t) = 0, khi đó P t at + C(t), a 6= 1. a. a = a−1 P cos a(t − 21 ) b. sin at = − + C(t), a 6= 2nπ. 2 sin a2 P sin a(t − 12 ) c. cos at = + C(t), a 6= 2nπ. a 2 sin 2 P d. log t = log Γ(t) + C(t), t > 0. P a ta+1 e. t = + C(t), a 6= −1. a+1 Chứng minh. a. Ta có ∆ at at+1 at + C(t) = − = at . a−1 a−1 a−1 14 Do đó X at a = + C(t), a 6= 1. a−1 t b. Ta có  cos a(t − 1 )   cos a(t − 1 )  2 2 ∆ − + C(t) = ∆ − a a 2 sin 2 2 sin 2 ∆ cos a(t − 12 ) =− 2 sin a2 cos a(t + 12 ) − cos a(t − 12 ) =− 2 sin a2 −2 sin a2 sin at =− 2 sin a2 = sin at. Như vậy X cos a(t − 12 ) sin at = − + C(t), a 6= 2nπ. 2 sin a2 c. Ta có  sin a(t − 1 )   sin a(t − 1 )  2 2 ∆ + C(t) = ∆ a a 2 sin 2 2 sin 2 ∆ sin a(t − 12 ) = 2 sin a2 sin a(t + 12 ) − sin a(t − 12 ) = 2 sin a2 2 sin a2 cos at = 2 sin a2 = cos at. Như vậy X sin a(t − 12 ) cos at = + C(t), a 6= 2nπ. 2 sin a2 15 d. Ta có     ∆ log Γ(t) + C(t) = ∆ log Γ(t) = log Γ(t + 1) − log Γ(t) = log Γ(t + 1) Γ(t) = log t. Như vậy X log t = log Γ(t) + C(t), t > 0. e. Ta có  ta+1  ta+1 ∆ + C(t) = ∆ a+1 a+1 (t + 1)a+1 − ta+1 = a+1 = log t. Như vậy X log t = log Γ(t) + C(t), t > 0. Ví dụ 1.4 Tìm nghiệm của phương trình sai phân sau y(t + 2) − 2y(t + 1) + y(t) = t(t − 1), t = 0, 1, 2, · · · thỏa mãn y(0) = −1, y(1) = 3. Phương trình sai phân trên có thể viết gọn dưới dạng ∆2 y(t) = t2 . Khi đó ta có t3 ∆y(t) = + C. 3 Và t4 y(t) = + Ct + D, 12 trong đó C, D là các hằng số. Sử dụng điều kiện y(0) = −1, y(1) = 3 ta tìm được D = −1 và C = 4. Do đó nghiệm duy nhất của phương trình là t4 y(t) = + 4t − 1. 12 16 Định lý 1.6 Tổng bất định có một số tính chất tổng quát sau đây P P P a. (y(t) + z(t)) = y(t) + z(t). P P b. (Dy(t)) = D y(t), trong đó D là hằng số. P P c. (y(t)∆z(t)) = y(t)z(t) − Ez(t)∆y(t). P P d. (Ey(t)∆z(t)) = y(t)z(t) − z(t)∆y(t). Chứng minh. a. b. Các ý a, b được suy trực tiếp từ tính tuyến tính của toán tử ∆. c. Ta có ∆(y(t)z(t)) = y(t)∆z(t) + Ez(t)∆y(t). Do vậy X  y(t)∆z(t) + Ez(t)∆y(t) = y(t)z(t) + C(t). Hay X (y(t)∆z(t)) = y(t)z(t) − X Ez(t)∆y(t). d. Ta có ∆(y(t)z(t) = ∆(z(t)y(t))) = z(t)∆y(t) + Ey(t)∆z(t). Do vậy  X z(t)∆y(t) + Ey(t)∆z(t) = y(t)z(t) + C(t). Hay X (Ey(t)∆z(t)) = y(t)z(t) − X z(t)∆y(t). Chú ý 1.2 Ta có một số chú ý sau đây: i. Tính chất c và d còn được biết đến với tên gọi "tổng từng phần." ii. Từ bây giờ về sau, ta sử dụng ký hiệu yn thay cho y(n) nếu tập xác định của y(n) là tập các số tự nhiên N = {0, 1, 2, 3, . . . }. 17 1.3. Biến đổi z Định nghĩa 1.6 Biến đổi z của dãy {yk } là một hàm Y (z) của biến phức được xác định bởi Y (z) = Z(yk ) = ∞ X yk k=0 zk , P∞ yk k=0 k z là hội tụ với |z| > R. Dãy {yk } được gọi là bị chặn mũ nếu tồn tại M > 0 và ta nói rằng biến đổi z là tồn tại nếu tồn tại số R > 0 sao cho và c > 1 sao cho |yk | ≤ M ck , ∀k ≥ 0. Định lý 1.7 Nếu dãy {yk } bị chặn mũ thì biến đổi z của {yk } là tồn tại. Chứng minh. Giả sử dãy {yk } bị chặn mũ, khi đó tồn tại M > 0 và c > 1 thỏa mãn |yk | ≤ M ck , ∀k ≥ 0. Ta có ∞ ∞ ∞ X X X |yk | c k yk ≤ M | | , | k| ≤ z |z|k z k=0 k=0 k=0 và tổng cuối hội tụ với |z| > c. Từ đó biến đổi z của dãy {yk } là tồn tại. Phần tiếp theo chúng tôi sẽ nêu lại một số định lý liên quan tới biến đổi z, chứng minh của các định lý này đã được trình bày rất chi tiết trong Chương 3, tài liệu [2]. Định lý 1.8 Nếu dãy {yk } bị chặn mũ thì nghiệm của phương trình sai phân cấp n là bị chặn mũ và do đó biến đổi z của nó là tồn tại. Định lý 1.9 (Tính tuyến tính) Nếu a và b là các hằng số thì Z(auk + bvk ) = aZuk + bZvk cho các z là xác định chung của U (z) và V (z). 18 Định lý 1.10 Nếu Y (z) = Z(yk ) với |z| > r thì n Z (k + n − 1) yk = n n nd Y (−1) z dz n (z), |z| > r. Định lý 1.11 Cho n là số nguyên dương n Z(yk+n ) = z Z(yk ) − n−1 X ym z n−m , m=0 Z(yk−n uk (n)) = z −n Z(yk ). Ví dụ 1.5 Tìm biến đổi z của dãy {yk = 1}. ∞ X 1 Y (z) = Z(1) = zk m=0 1 1 − z −1 z = , |z| > 1. z−1 = Ví dụ 1.6 Tìm Z(k 2 ). Z(k 2 ) = Z(k.k) d = −z Z(k) dz  d z = −z dz (z − 1)2 z(z + 1) = . (z − 1)3
- Xem thêm -

Tài liệu liên quan

Tài liệu vừa đăng

Tài liệu xem nhiều nhất