Tài liệu Phương pháp cim đối với bài toán biên elliptic có bước nhảy gián đoạn qua mặt phân cách

ĐẠI HỌC THÁI NGUYÊN TRƯỜNG ĐẠI HỌC CÔNG NGHỆ THÔNG TIN VÀ TRUYỀN THÔNG ĐÀO THỊ THẮM PHƯƠNG PHÁP CIM ĐỐI VỚI BÀI TOÁN BIÊN ELLIPTIC CÓ BƯỚC NHẢY GIÁN ĐOẠN QUA MẶT PHÂN CÁCH LUẬN VĂN THẠC SĨ KHOA HỌC MÁY TÍNH Thái Nguyên – 2011 ĐẠI HỌC THÁI NGUYÊN Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn TRƯỜNG ĐẠI HỌC CÔNG NGHỆ THÔNG TIN VÀ TRUYỀN THÔNG ĐÀO THỊ THẮM PHƯƠNG PHÁP CIM ĐỐI VỚI BÀI TOÁN BIÊN ELLIPTIC CÓ BƯỚC NHẢY GIÁN ĐOẠN QUA MẶT PHÂN CÁCH Chuyên ngành : Khoa học Máy tính Mã số : 60 48 01 LUẬN VĂN THẠC SĨ KHOA HỌC MÁY TÍNH NGƯỜI HƯỚNG DẪN KHOA HỌC: TS. VŨ VINH QUANG Thái Nguyên – năm 2011 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn LỜI CAM ĐOAN Tôi xin cam đoan Luận văn “ Phương pháp CIM đối với bài toán biên elliptic có bước nhảy gián đoạn qua mặt phân cách” là công trình nghiên cứu của riêng tôi dưới sự hướ ng dẫn của TS . Vũ Vinh Quang . Tôi xin chị u trách nhiệm về lời cam đoan của mì nh. Các số liệu và thông tin sử dụng trong luận văn này là trung thực. Tác giả Đào Thị Thắm Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn MỤC LỤC LỜI CAM ĐOAN MỤC LỤC MỞ ĐẦU ....................................................................................................................1 Chương 1 ....................................................................................................................3 Các kiến thức cơ bản.................................................................................................3 1.1. Các khái niệm về phương trình đạo hàm riêng .............................................3 1.1.1. Khái niệm phương trình đạo hàm riêng ..................................................3 1.1.2. Một số phương trình đạo hàm riêng tiêu biểu. .......................................4 1.1.3. Phân loại phương trình cấp hai tuyến tính ..............................................5 1.2. Phương pháp lưới giải phương trình đạo hàm riêng......................................8 1.2.1. Bài toán vi phân ......................................................................................8 1.2.2. Lưới sai phân .......................................................................................... 8 1.2.3. Hàm lưới .................................................................................................9 1.2.4. Đạo hàm lưới .......................................................................................... 9 1.2.5. Bài toán sai phân ...................................................................................10 1.2.6. Bài toán biên elliptic .............................................................................10 Chương 2 ..................................................................................................................28 Phương pháp CIM (Coupling Interface Method) ................................................28 2.1. Giới thiệu về bài toán biên với mặt phân cách gián đoạn.............................. 28 2.2. Cơ sở các phương pháp CIM đối với phương trình elipptic ........................... 29 2.2.1. Trường hợp không gian một chiều ........................................................... 30 2.2.2. 2.3. Trường hợp không gian hai chiều ......................................................... 38 Một số kết quả thực nghiệm với CIM2 .......................................................45 Chương 3 ..................................................................................................................49 Hướng tiếp cận phương pháp chia miền đối với bài toán biên gián đoạn qua mặt phân cách ..........................................................................................................49 3.1. Cơ sở phương pháp chia miền .....................................................................49 3.2. Bài toán biên với hệ số đạo hàm gián đoạn .................................................52 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 3.3. Bài toán biên gián đoạn qua mặt phân cách ................................................54 3.4. Các kết quả thực nghiệm .............................................................................56 KẾT LUẬN ..............................................................................................................65 DANH MỤC TÀI LIỆU THAM KHẢO ............................................................... 66 PHỤ LỤC .................................................................................................................68 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn MỞ ĐẦU Trong trường hợp tổng quát, việc tìm nghiệm đúng của lớp các bài toán biên mà chủ yếu là phương trình elliptic cấp hai là không thực hiện được. Việc nghiên cứu giải gần đúng các bài toán biên mà tiêu biểu là các bài toán được biểu diễn bằng các phương trình cấp hai đã và đang là một lĩnh vực rất quan tâm của các nhà toán học. Trong nhiều năm qua đã có nhiều công trình nghiên cứu về lĩnh vực này, mục đích chính của các phương pháp là đưa bài toán vi phân về bài toán rời rạc trên một điểm lưới. Tuy nhiên khi miền hình học là miền phức tạp, dữ liệu hoặc các hệ số của phương trình là gián đoạn thì việc áp dụng một phương pháp số nào đó cho cả miền gặp nhiều khó khăn. Để giải quyết khó khăn trên, trong nhiều năm qua đã có nhiều công trình nghiên cứu về lĩnh vực này. Các hướng nghiên cứu chủ yếu là đưa ra các phương pháp sai phân đặc biệt xung quanh lân cận các điểm kỳ dị hoặc biên phân chia để đưa bài toán đang xét về các hệ phương trình sai phân và việc tìm nghiệm bằng số của bài toán chuyển về việc giải các hệ phương trình đại số bằng các phương pháp đúng hoặc gần đúng. Một hướng thứ hai là sử dụng phương pháp chia miền chuyển bài toán trên miền đang xét về hai bài toán không chứa các điểm kỳ dị, sau đó xuất phát từ lời giải các bài toán trên hai miền ta thu được nghiệm của bài toán gốc. Nội dung chính của luận văn là đặt vấn đề tìm hiểu phương pháp sai phân đặc biệt được gọi là phương pháp CIM (Coupling interface method) giải phương trình elliptic cấp hai trong trường hợp miền đang xét tồn tại mặt phân cách tại đó xảy ra gián đoạn của hàm và đạo hàm và đồng thời nghiên cứu phương pháp chia miền giải bài toán trên, tiến hành tính toán thử nghiệm và so sánh giữa phương pháp chia miền và phương pháp CIM. Luận văn cấu trúc gồm 3 chương: Chương 1: Đưa ra các kiến thức cơ bản về phương trình đạo hàm riêng, cơ sở phương pháp lưới, thuật toán thu gọn khối lượng tính toán giải hệ phương trình lưới Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 1 và giới thiệu thư viện chương trình giải phương trình elliptic với hệ số hằng số trong miền chữ nhật. Chương 2: Trình bày cơ sở phương pháp CIM bao gồm: phương pháp CIM1, CIM2 trong không gian một chiều và không gian hai chiều xây dựng phương pháp sai phân đặc biệt đối với bài toán biên có mặt gián đoạn, các thuật toán cơ bản về các phương pháp tương ứng, các kết quả thực nghiệm đối với các bài toán cụ thể. Chương 3: Trình bày cơ sở phương pháp chia miền, xây dựng các sơ đồ lặp giải bài toán biên elliptic tồn tại mặt gián đoạn theo hướng hiệu chỉnh giá trị hàm trên biên, xây dựng các chương trình thực nghiệm trong các bài toán cụ thể và bài toán tổng quát. So sánh phương pháp chia miền với phương pháp CIM. Phần cuối của luận văn là các chương trình thực nghiệm được viết trên ngôn ngữ Matlab version7.0. Tác giả xin bày tỏ lòng biết ơn sâu sắc đối với thầy giáo hướng dẫn TS. Vũ Vinh Quang đã tận tình hướng dẫn tác giả trong suốt quá trình làm luận văn. Tác giả xin chân thành cảm ơn các Thầy, Cô giáo Viện Công nghệ thông tin, Trường Đại học Công nghệ thông tin và Truyền thông – Đại học Thái Nguyên đã tham gia giảng dạy, giúp đỡ tác giả trong suốt quá trình học tập và nghiên cứu. Mặc dù đã cố gắng rất nhiều song nội dung bản luận văn không tránh khỏi những thiếu sót. Tác giả rất mong nhận được sự chỉ bảo, đóng góp ý kiến của các Thầy Cô và các bạn để luận văn thêm hoàn thiện. Tác giả Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 2 Chương 1 Các kiến thức cơ bản Trong Chương này luận văn trình bày các kiến thức cơ bản bao gồm: các khái niệm cơ bản về phương trình đạo hàm riêng, phương pháp lưới giải phương trình đạo hàm riêng, thuật toán thu gọn khối lượng tính toán giải hệ phương trình vectơ ba điểm và đặc biệt là giới thiệu thư viện RC2009 giải số bài toán biên elliptic với hệ số hằng. Các kiến thức cơ bản được tham khảo trong các tài liệu [1, 2, 3, 8, 10]. 1.1. Các khái niệm về phương trình đạo hàm riêng 1.1.1. Khái niệm phương trình đạo hàm riêng Phương trình đạo hàm riêng là một lĩnh vực quan trọng của toán học. Có rất nhiều mô hình trong tự nhiên được mô tả bởi một phương trình hay một hệ phương trình vi phân nói chung và phương trình đạo hàm riêng nói riêng. Với hàm số một biến 𝑦 = 𝑦(𝑥) ta có khái niệm đạo hàm 𝑦 ′ (𝑥): 𝑦 𝑥 + ∆𝑥 − 𝑦(𝑥) ∆𝑥→0 ∆𝑥 𝑦 ′ 𝑥 = lim khái niệm phương trình vi phân 𝑦 ′ = 𝑓 𝑥, 𝑦 và khái niệm bài toán Cauchy: Tìm hàm số 𝑦 = 𝑦(𝑥) xác định tại 𝑥 ∈ [𝑥0 , 𝑋] sao cho: 𝑦 ′ = 𝑓 𝑥, 𝑦 , 𝑥0 < 𝑥 ≤ 𝑋, 𝑦 𝑥0 = 𝜂 Trong đó: 𝑓 𝑥, 𝑦 là hàm số cho trước. 𝑥0 , 𝑋, 𝜂 là những số cho trước. Với hàm số nhiều biến số ta cũng gặp những khái niệm và bài toán tương tự. Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 3 Xét bài toán hai biến số 𝑢 = 𝑢(𝑥, 𝑦), ta có đạo hàm riêng cấp 1 đối với biến 𝑥: 𝜕𝑢 𝑢 𝑥 + ∆𝑥, 𝑦 − 𝑢(𝑥, 𝑦) = lim 𝜕𝑥 ∆𝑥→0 ∆𝑥 đạo hàm riêng cấp 1 đối với biến y: 𝜕𝑢 𝑢 𝑥, 𝑦 + ∆𝑦 − 𝑢(𝑥, 𝑦) = lim 𝜕𝑦 ∆𝑦→0 ∆𝑦 và các đạo hàm riêng cấp 2: 𝜕2 𝑢 𝜕 𝜕𝑢 = , 𝜕𝑥 2 𝜕𝑥 𝜕𝑥 𝜕2 𝑢 𝜕 𝜕𝑢 = , 𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕2 𝑢 𝜕 𝜕𝑢 = 𝜕𝑦 2 𝜕𝑦 𝜕𝑦 𝜕2 𝑢 𝜕 𝜕𝑢 = 𝜕𝑦𝜕𝑥 𝜕𝑥 𝜕𝑦 Nếu các đạo hàm riêng 𝜕2𝑢 𝜕𝑥𝜕𝑦 và 𝜕2𝑢 𝜕𝑦𝜕𝑥 là các hàm liên tục thì chúng bằng nhau. Phương trình: 𝜕2 𝑢 𝜕2 𝑢 𝜕2 𝑢 𝜕𝑢 𝜕𝑢 𝐴 𝑥, 𝑦 + 𝐵 𝑥, 𝑦 + 𝐶 𝑥, 𝑦 + 𝐷 𝑥, 𝑦 + 𝐸 𝑥, 𝑦 + 2 2 𝜕𝑥 𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝐹 𝑥, 𝑦 𝑢 = 𝑓(𝑥, 𝑦) là phương trình đạo hàm riêng cấp 2. 1.1.2. Một số phương trình đạo hàm riêng tiêu biểu. 1. Phương trình Laplace do Laplace đưa ra vào năm 1780 𝑛 𝑢𝑥 𝑖 𝑥 𝑖 = 0, 𝑥 ∈ 𝑅𝑛 ∆𝑢 = 𝑖=1 2. Phương trình Helmholtz được Helmholtz nghiên cứu vào năm 1860 −∆𝑢 = 𝜆𝑢 3. Phương trình chuyển dịch tuyến tính Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 4 𝑛 𝑢𝑡 + 𝑏𝑖 𝑢𝑥 𝑖 = 0 𝑖=1 4. Phương trình Liouville được nghiên cứu vào khoảng năm 1851 𝑛 𝑢𝑡 − (𝑏𝑖 𝑢)𝑥 𝑖 = 0 𝑖=1 5. Phương trình truyền nhiệt được Fourier công bố vào năm 1810 – 1822 𝑢𝑡 = ∆𝑢 6. Phương trình Schrodinger (1926) 𝑖𝑢𝑡 + ∆𝑢 = 0 7. Phương trình truyền sóng được D’Alembert đưa ra năm 1752 𝑢𝑡𝑡 − ∆𝑢 = 0 và dạng tổng quát của nó là: 𝑛 𝑢𝑡𝑡 − 𝑛 𝑎𝑖𝑗 𝑢𝑥 𝑖 𝑥 𝑗 + 𝑖=1 𝑏𝑖 𝑢𝑥 𝑖 = 0 𝑖=1 1.1.3. Phân loại phương trình cấp hai tuyến tính Giả sử 𝑢 = 𝑢(𝑝, 𝑞) là hàm số của hai biến độc lập 𝑝, 𝑞. Kí hiệu: 𝜕𝑢 𝜕𝑢 𝜕2 𝑢 𝜕2 𝑢 𝜕2 𝑢 𝑢𝑝 = ; 𝑢𝑞 = ; 𝑢𝑝𝑝 = 2 ; 𝑢𝑞𝑞 = 2 ; 𝑢𝑝𝑞 = 𝑢𝑞𝑝 = ; 𝜕𝑝 𝜕𝑞 𝜕𝑝 𝜕𝑞 𝜕𝑝𝜕𝑞 Xét phương trình đạo hàm riêng cấp 2 á tuyến: 𝐴𝑢𝑝𝑝 + 2𝐵𝑢𝑝𝑞 + 𝐶𝑢𝑞𝑞 = 𝐹 (1.1) Với 𝐴, 𝐵, 𝐶, 𝐹 là những hàm số phụ thuộc 𝑝, 𝑞, 𝑢𝑝 , 𝑢𝑞 . Giả sử phương trình trên có nghiệm là 𝑢 = 𝑢(𝑝, 𝑞) đủ trơn. Xét 𝛤 là một đường cong nào đó của mặt phẳng (𝑝, 𝑞) nằm trong miền xác định của hàm 𝑢(𝑝, 𝑞) và có phương trình 𝑞 = 𝑞(𝑝), hay 𝜑 𝑝, 𝑞 = 0. Ta có: 𝑑 𝑢𝑝 = 𝑢𝑝𝑝 𝑑𝑝 + 𝑢𝑝𝑞 𝑑𝑞 ; 𝑑 𝑢𝑞 = 𝑢𝑞𝑝 𝑑𝑝 + 𝑢𝑞𝑞 𝑑𝑞 Vậy ta có hệ: Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 5 𝐴𝑢𝑝𝑝 + 2𝐵𝑢𝑝𝑞 + 𝐶𝑢𝑞𝑞 = 𝐹(𝑝, 𝑞, 𝑢, 𝑢𝑝 , 𝑢𝑞 ) 𝑢𝑝𝑝 𝑑𝑝 + 𝑢𝑝𝑞 𝑑𝑞 = 𝑑 𝑢𝑝 𝑢𝑝𝑝 𝑑𝑝 + 𝑢𝑝𝑞 𝑑𝑞 = 𝑑 𝑢𝑝 hay ở dạng ma trận: 𝐴 2𝐵 𝐶 𝑑𝑝 𝑑𝑞 0 0 𝑑𝑝 𝑑𝑞 𝑢𝑝𝑝 𝑢𝑝𝑞 𝑢𝑞𝑞 𝑓 = 𝑑 𝑢𝑝 𝑑(𝑢𝑞 ) Hệ này luôn có nghiệm vì ta đã giả sử phương trình (1.1) có nghiệm 𝑢 = 𝑢(𝑝, 𝑞) đủ trơn. Xét ma trận của hệ: 𝐴 2𝐵 𝐶 𝑀 = 𝑑𝑝 𝑑𝑞 0 0 𝑑𝑝 𝑑𝑞 Nếu 𝐷𝑒𝑡(𝑀) ≠ 0 trên 𝛤 thì hệ trên có nghiệm duy nhất trên 𝛤, nghĩa là trên 𝛤 các đạo hàm cấp hai của 𝑢 được xác định một cách duy nhất theo vế phải. Nếu 𝐷𝑒𝑡 𝑀 = 0 trên 𝛤 thì hệ trên vẫn có nghiệm duy nhất trên 𝛤 vì ta đã xuất phát từ giả thiết phương trình (1.1) có nghiệm 𝑢, nhưng nghiệm đó không duy nhất nữa, nghĩa là trên 𝛤 các đạo hàm cấp hai của 𝑢 xác định một cách không duy nhất theo vế phải. Trong trường hợp này ta gọi 𝛤 là một “đường đặc trưng” của phương trình đạo hàm riêng (1.1). Như vậy đường đặc trưng xác định bởi điều kiện 𝐷𝑒𝑡 𝑀 = 0, điều kiện này viết như sau: 𝐷𝑒𝑡 𝑀 = 𝐴(𝑑𝑞)2 − 2𝐵𝑑𝑝𝑑𝑞 + 𝐶(𝑑𝑝)2 = 0 (1.2) hay: 𝐴( 𝑑𝑞 2 𝑑𝑞 ) − 2𝐵 +𝐶 =0 𝑑𝑝 𝑑𝑝 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên (1.3) http://www.lrc-tnu.edu.vn 6 Trong đó 𝑑𝑞 𝑑𝑝 𝑑𝑞 𝑑𝑝 là hệ số góc của tiếp tuyến của đường đặc trưng, người ta gọi là phương đặc trưng tại điểm (𝑝, 𝑞). Vậy phương trình (1.3) xác định phương đặc trưng, nó là phương trình vi phân của đường đặc trưng. Phương trình 1.3 là một phương trình bậc hai đối với 𝑑𝑞 𝑑𝑝 Xét ∆= 𝐵2 − 𝐴𝐶  Nếu 𝐵2 − 𝐴𝐶 > 0 tại (𝑝, 𝑞) ∈ miền 𝛺 nào đó thì phương trình (1.3) có hai nghiệm thực khác nhau tại (𝑝, 𝑞) ∈ 𝛺: 𝑑𝑞 𝐵 ± −𝐵2 + 𝐴𝐶 = 𝑑𝑝 𝐴 Khi đó tại mỗi (𝑝, 𝑞) ∈ 𝛺 có hai phương đặc trưng thực khác nhau, ta nói phương trình (1.1) thuộc loại hypebol trong 𝛺.  Nếu 𝐵2 − 𝐴𝐶 = 0 tại (𝑝, 𝑞) ∈ miền 𝛺 nào đó thì phương trình (1.3) có hai nghiệm thực trùng nhau tại (𝑝, 𝑞) ∈ 𝛺: 𝑑𝑞 𝐵 = 𝑑𝑝 𝐴 Khi đó tại mỗi (𝑝, 𝑞) ∈ 𝛺 có hai phương đặc trưng thực trùng nhau, ta nói phương trình (1.1) thuộc loại parabol trong 𝛺.  Nếu 𝐵2 − 𝐴𝐶 < 0 tại (𝑝, 𝑞) ∈ miền 𝛺 nào đó thì phương trình (1.3) không có nghiệm thực nào mà chỉ có hai nghiệm phức liên hợp tại (𝑝, 𝑞) ∈ 𝛺: 𝑑𝑞 𝐵 ± 𝑖 −𝐵2 + 𝐴𝐶 = 𝑑𝑝 𝐴 Khi đó tại mỗi (𝑝, 𝑞) ∈ 𝛺 không có phương đặc trưng thực nào mà chỉ có hai phương đặc trưng ảo liên hợp, ta nói phương trình (1.1) thuộc loại elip trong 𝛺. Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 7 Ví dụ: 𝜕2 𝑢 𝜕2 𝑢 Phương trình Laplace: ∆𝑢 = 2 + 2 = 0 𝜕𝑥 𝜕𝑦 𝜕2 𝑢 𝜕2 𝑢 và phương trình Poisson: ∆𝑢 = 2 + 2 = 𝑄(𝑥, 𝑦) 𝜕𝑥 𝜕𝑦 có 𝐴 = 𝐶 = 1, 𝐵 = 0, do đó 𝐵2 − 𝐴𝐶 = 0 − 1.1 < 0 tại mọi (𝑥, 𝑦), nên là các phương trình elip trong cả hai mặt phẳng (𝑥, 𝑦). 1.2. Phương pháp lưới giải phương trình đạo hàm riêng Phương pháp lưới (hay còn gọi là phương pháp sai phân) là phương pháp được áp dụng rộng rãi trên nhiều lĩnh vực khoa học, kỹ thuật. Nội dung chính của nó là đưa bài toán vi phân đang xét về việc giải hệ phương trình sai phân (tức là hệ thức hoặc các hệ thức liên hệ các giá trị của các hàm số tại các điểm khác nhau) bằng các phương pháp đại số. 1.2.1. Bài toán vi phân Cho hai số 𝑎 và 𝑏 với 𝑎 < 𝑏. Tìm hàm 𝑢 = 𝑢(𝑥) xác định tại 𝑎 < 𝑥 < 𝑏 thỏa mãn: 𝐿𝑢 = − 𝑘𝑢 ′ ′ 𝑢 𝑎 = 𝛼, 𝑢 𝑏 =𝛽 + 𝑞𝑢 = 𝑓 𝑥 với 𝑎 < 𝑥 < 𝑏 (1.4) (1.5) trong đó 𝑘 = 𝑘 𝑥 , 𝑞 = 𝑞 𝑥 , 𝑓(𝑥) là những hàm số cho trước đủ trơn thỏa mãn: 0 < 𝑐0 ≤ 𝑘 𝑥 ≤ 𝑐1 , 𝑐0 , 𝑐1 = 𝑐𝑜𝑛𝑠𝑡, 𝑞(𝑥) ≥ 0 (1.6) 𝛼, 𝛽 là những số cho trước. Giả sử bài toán (1.4) – (1.5) có nghiệm duy nhất 𝑢 đủ trơn trên [𝑎, 𝑏] Đây chính là bài toán biên của phương trình elip một chiều. 1.2.2. Lưới sai phân Xét bài toán Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 8 −∆𝑢 = 𝑓, 𝑥 ∈ 𝛺, 𝑢 = 𝑔, 𝑥 ∈ 𝜕𝛺 (1.7) trong đó 𝛺 = {(𝑥, 𝑦) ∈ 𝑅2 , 𝑎 ≤ 𝑥 ≤ 𝑏, 𝑐 ≤ 𝑦 ≤ 𝑑}, chọn 2 số nguyên 𝑁 > 1 và 𝑀 > 1, đặt 𝑕 = (𝑏 − 𝑎)/𝑁 gọi là bước lưới theo 𝑥, 𝑘 = 𝑑 − 𝑐 /𝑀 gọi là bước lưới theo 𝑦. Đặt 𝑥𝑖 = 𝑎 + 𝑖𝑕, 𝑦𝑗 = 𝑐 + 𝑗𝑕, 𝑖 = 0. . 𝑁, 𝑗 = 0. . 𝑀. Mỗi điểm (𝑥𝑖 , 𝑦𝑗 ) gọi là một nút lưới ký hiệu là nút (𝑖, 𝑗). Tập tất cả các nút trong ký hiệu là 𝛺𝑕𝑘 . Nút ở trên biên 𝛤 gọi là nút biên; tập tất cả các nút biên ký hiệu là 𝛤𝑕𝑘 , tập 𝛺𝑕𝑘 = 𝛺𝑕𝑘 ∪ 𝛤𝑕𝑘 gọi là một lưới sai phân trên 𝛺−. 1.2.3. Hàm lưới Những hàm số xác định tại các nút của lưới 𝛺𝑕 được gọi là hàm lưới. Giá trị của hàm lưới 𝑣 tại nút 𝑥𝑖 viết là 𝑣𝑖 . Một hàm số 𝑢(𝑥) xác định tại mọi 𝑥 ∈ [𝑎, 𝑏] sẽ tạo ra hàm lưới 𝑢 có giá trị tại nút 𝑥𝑖 là 𝑢𝑖 = 𝑢(𝑥𝑖 ). 1.2.4. Đạo hàm lưới Xét hàm lưới 𝑣. Đạo hàm lưới tiến cấp một của v, ký hiệu là 𝑣𝑥 , có giá trị tại nút 𝑥𝑖 là: 𝑣𝑥 𝑖 = 𝑣𝑖+1 −𝑣𝑖 𝑕 Đạo hàm lưới lùi cấp một của 𝑣, ký hiệu là 𝑣𝑥 , có giá trị tại nút 𝑥𝑖 là: 𝑣𝑥 𝑖 = 𝑣𝑖 −𝑣𝑖−1 𝑕 Khi 𝑕 bé thì đạo hàm lưới “xấp xỉ” được đạo hàm thường. Đạo hàm lưới cấp hai 𝑣𝑥 𝑥 : 𝑣𝑥 𝑥 = 𝑣𝑥 𝑖 +1 − 𝑣𝑥 𝑖 𝑣𝑖+1 −𝑣𝑖 𝑣𝑖 −𝑣𝑖−1 𝑣𝑖+1 −2𝑣𝑖 + 𝑣𝑖−1 = − /𝑕 = 𝑕 𝑕 𝑕 𝑕2 Nếu 𝑎 là một hàm lưới nữa thì: Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 9 𝑎𝑣𝑥 𝑥 = 𝑎𝑖+1 𝑣𝑥 𝑖 +1 − 𝑎𝑖 𝑣𝑥 𝑖 𝑎𝑖+1 𝑣𝑖+1 −(𝑎𝑖+1 + 𝑎𝑖 )𝑣𝑖 + 𝑎𝑖 𝑣𝑖−1 = 𝑕 𝑕2 1.2.5. Bài toán sai phân Ta tìm cách tính gần đúng giá trị của nghiệm đúng 𝑢(𝑥𝑖 ) tại các nút 𝑥𝑖 ∈ 𝛺𝑕 . Gọi các giá trị gần đúng đó là 𝑣𝑖 . Để tìm 𝑣𝑖 ta thay bài toán (1.4) – (1.5) bởi bài toán sai phân: 𝐿𝑕 𝑣 ≡ − 𝑎𝑣𝑥 𝑥𝑖 + 𝑞𝑖 𝑣𝑖 = 𝑓𝑖 𝑣0 = 𝛼, 𝑣𝑁 = 𝛽 Trong đó: 𝑎𝑖 = 𝑘 𝑥𝑖 − 𝑕 , 𝑞𝑖 = 𝑞 𝑥𝑖 , 2 𝑓𝑖 = 𝑓 𝑥𝑖 1.2.6. Bài toán biên elliptic Xét bài toán biên: 𝜕2 𝑢 𝜕2 𝑢 + 𝑘 𝑥 + 𝑐 𝑥 𝑢 = 𝑓 𝑥 ,𝑥 ∈ 𝛺 2 𝜕𝑥12 𝜕𝑥22 𝑙𝑢 = 𝑔 𝑥 , 𝑥 ∈ 𝜕𝛺 𝑘1 𝑥 (1.8) trong đó: 𝛺 là miền hình chữ nhật có kích thước hai cạnh là 𝐿1 và 𝐿2 , 𝑙 là toán tử điều kiện biên, 𝑓(𝑥) và 𝑔(𝑥) là các hàm số cho trước. Xét trường hợp tổng quát, với điều kiện biên 𝑙𝑢 = 𝑔 𝑥 là điều kiện biên Dirichlet hoặc Neumann trên các phần biên khác nhau trong đó tồn tại ít nhất một điều kiện biên Dirichlet trên một cạnh để đảm bảo bài toán có nghiệm duy nhất. Trong luận văn sẽ sử dụng phương pháp sai phân để nghiên cứu việc giải số bài toán elliptic cấp hai đã cho ở trên. Đưa vào không gian lưới 𝛺𝑘𝑕 = 𝑥𝑖𝑗 = 𝑖𝑘, 𝑗𝑕 , 𝑖 = 0, 𝑀 , 𝑗 = 0, 𝑁 với 𝑘 = 𝐿1 𝐿2 ,𝑕 = 𝑀 𝑁 Khi đó bài toán vi phân đang xét luôn được đưa về các hệ phương trình véctơ ba điểm có dạng Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 10 −𝑌𝑗 −1 + 𝐶𝑌𝑗 − 𝑌𝑗 +1 , 𝑗 = 1, 𝑁 − 1, 𝑌0 = 𝐹0 , 𝑌𝑁 = 𝐹𝑁 (1.9) trong trường hợp bài toán biên Dirichlet và có dạng 𝐶𝑌0 − 2𝑌1 = 𝐹0 −𝑌𝑗 −1 + 𝐶𝑌𝑗 − 𝑌𝑗 +1 , 𝑗 = 1, 𝑁 − 1, −2𝑌𝑁−1 + 𝐶𝑌𝑁 = 𝐹𝑁 (1.10) trong trường hợp bài toán Neumann, trong đó 𝐶 là ma trận ba đường chéo trội, 𝑌𝑗 là các véc tơ nghiệm, các véc tơ 𝐹𝑗 chứa các giá trị hàm vế phải và giá trị hàm hoặc đạo hàm trên biên. Như vậy, để giải được bài toán (1.8) bằng phương pháp số, điều quan trọng nhất là ta phải xác định được thuật toán nhanh giải các hệ phương trình véc tơ ba điểm (1.9), (1.10) là các hệ phương trình đại số tuyến tính. Có nhiều phương pháp khác nhau để giải được các hệ trên. Tuy nhiên do tính chất đặc biệt của hệ, trong luận văn tác giả giới thiệu phương pháp thu gọn khối lượng tính toán của SamarskijNicolaev đề xuất [4] với độ phức tạp tính toán 𝑂(𝑀𝑁𝑙𝑜𝑔𝑁) và các kết quả xây dựng thư viện chương trình tìm nghiệm bằng số của bài toán elliptic với các điều kiện biên khác nhau được thiết kế trên môi trường Matlab. 1.2.6.1. Phương pháp thu gọn khối lượng tính toán a. Bài toán biên thứ nhất: xét phương trình (1.9) Ý tưởng của phương pháp: Thực hiện khử liên tiếp các ẩn 𝑌𝑗 (với các 𝑗 lẻ), tiếp theo từ các phương trình còn lại khử các 𝑌𝑗 với 𝑗 là bội của 2, bội của 4,… Mỗi bước khử sẽ giảm được một nửa số ẩn. Như vậy, nếu 𝑁 = 2𝑛 thì sau một số lần khử sẽ chỉ còn lại một phương trình chứa véc tơ ẩn 𝑌𝑁/2 mà từ đó 𝑌𝑁/2 có thể tính được qua 𝑌0 và 𝑌𝑁 . Sau khi đã có được 𝑌0 , 𝑌𝑁/2 và 𝑌𝑁 thì quá trình ngược lại là tìm các 𝑌𝑗 với j là bội của 𝑁/4, bội của 𝑁/8,… Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 11 Cụ thể như sau: Giả sử 𝑁 = 2𝑛 , 𝑛 > 0. Ký hiệu 𝐶 (0) = 𝐶, 𝐹𝑗 (0) = 𝐹𝑗 , 𝑗 = 0, 𝑁 − 1 Khi đó (1.9) được viết dưới dạng: −𝑌𝑗 −1 + 𝐶 (0) 𝑌𝑗 − 𝑌𝑗 +1 = 𝐹𝑗 , 𝑗 = 0, 𝑁 − 1 𝑌0 = 𝐹0 𝑌𝑁 = 𝐹𝑁 Bước khử thứ nhất: Từ các phương trình của (1.9), khử các 𝑌𝑗 với 𝑗 lẻ. Muốn vậy ta viết 3 phương trình liên tiếp trong (1.9). (0) −𝑌𝑗 −2 + 𝐶 (0) 𝑌𝑗 −1 − 𝑌𝑗 = 𝐹𝑗 −1 ; −𝑌𝑗 −1 + 𝐶 0 𝑌𝑗 − 𝑌𝑗 +1 = 𝐹𝑗 0 ; (0) −𝑌𝑗 + 𝐶 (0) 𝑌𝑗 +1 − 𝑌𝑗 +2 = 𝐹𝑗 +1 . Nhân hai vế của phương trình thứ 2 với 𝐶 0 , sau đó cộng cả 3 phương trình lại ta được: −𝑌𝑗 −2 + 𝐶 𝑌0 = 𝐹0 𝑌𝑁 = 𝐹𝑁 1 1 𝑌𝑗 − 𝑌𝑗 +2 = 𝐹𝑗 −1 ; 𝑗 = 2,4, … , 𝑁 − 2 (1.11) trong đó 𝐶 (1) = (𝐶 (0) )2 − 2𝐸; 𝐹𝑗 (1) (0) = 𝐹𝑗 −1 + 𝐶 (0) 𝐹𝑗 (0) (0) + 𝐹𝑗 +1 ; 𝑗 = 2,4, … , 𝑁 − 2 Hệ (1.11) chỉ chứa các 𝑌𝑗 với j chẵn, số véc tơ ẩn 𝑌𝑗 là 𝑁 2 −1. Nếu giải được hệ này thì các 𝑌𝑗 với 𝑗 lẻ sẽ tìm được từ phương trình: 𝐶 (0) 𝑌𝑗 = 𝐹𝑗 (0) + 𝑌𝑗 −1 + 𝑌𝑗 +1 , 𝑗 = 1, 3, … , 𝑁 − 1 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên (1.12) http://www.lrc-tnu.edu.vn 12 Bước khử thứ hai: Tiến hành khử các 𝑌𝑗 của hệ (1.11) với 𝑗 là bội của 2 nhưng không là bội của 4. Muốn vậy ta viết 3 phương trình liên tiếp của (1.11). 1 −𝑌𝑗 −4 + 𝐶 1 𝑌𝑗 −2 − 𝑌𝑗 = 𝐹𝑗 −2 −𝑌𝑗 −2 + 𝐶 1 𝑌𝑗 − 𝑌𝑗 +2 = 𝐹𝑗 −𝑌𝑗 + 𝐶 1 1 1 𝑌𝑗 +2 − 𝑌𝑗 +4 = 𝐹𝑗 +2 𝑗 = 4, 8, … , 𝑁 − 4 Nhân hai vế của phương trình thứ hai với 𝐶 1 vào bên trái rồi cộng cả 3 phương trình ta thu được: −𝑌𝑗 −4 + 𝐶 𝑌0 = 𝐹0 𝑌𝑁 = 𝐹𝑁 𝑌𝑗 − 𝑌𝑗 +4 = 𝐹𝑗 2 , 𝑗 = 4, 8, … , 𝑁 − 4 2 (1.13) trong đó: 𝐶 2 = 𝐶 1 2 − 2𝐸 1 𝐹𝑗 2 = 𝐹𝑗 −2 +𝐶 hệ (1.13) chỉ chứa 𝑁 4 1 1 𝐹𝑗 1 + 𝐹𝑗 +2 , 𝑗 = 4, 8, … , 𝑁 − 4 −1 véc tơ ẩn 𝑌𝑗 , trong đó j là bội của 4. Nếu giải được hệ này thì các 𝑌𝑗 với j là bội của 2 nhưng không là bội của 4 sẽ tìm được từ phương trình: 𝐶 1 𝑌𝑗 = 𝑌𝑗 −2 + 𝑌𝑗 +2 + 𝐹𝑗 1 , 𝑗 = 2, 6, 10, … , 𝑁 − 2. (1.14) Tiếp tục quá trình khử này. Kết quả là sau bước khử thứ l ta nhận được hệ gồm 𝑁 2𝑙 −1 ẩn 𝑌𝑗 , trong đó 𝑗 là bội của 2𝑙 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 13 𝑆ố 𝑕ó𝑎 𝑏ở𝑖 𝑇𝑟𝑢𝑛𝑔 𝑡â𝑚 𝐻ọ𝑐 𝑙𝑖ệ𝑢 – Đạ𝑖 𝑕ọ𝑐 𝑇𝑕á𝑖 𝑁𝑔𝑢𝑦ê𝑛 1 1 𝑡𝑛𝑢. 𝑒𝑑𝑢. 𝑣𝑛 𝑕𝑡𝑡𝑝://𝑤𝑤𝑤. 𝑙𝑟𝑐 − 𝑙 𝑙 𝑙 −𝑌𝑗−2𝑙 + 𝐶 𝑌𝑗 − 𝑌𝑗+2𝑙 = 𝐹𝑗 , 𝑗 = 2 , 2. 2 , 3. 2 , … , 𝑁 − 2 𝑙 𝑌0 = 𝐹0 𝑌𝑁 = 𝐹𝑁 (1.15) và nhóm các phương trình 𝐶 𝑘−1 𝑌𝑗 = 𝐹𝑗 𝑘−1 + 𝑌𝑗 −2𝑘−1 + 𝑌𝑗 +2𝑘−1 , 𝑗 = 2𝑘−1 , 3. 2𝑘−1 , 5. 2𝑘−1 , … , 𝑁 − 2𝑘−1 , 𝑘 = 𝑙, 𝑙 − 1, … , 1 trong đó các ma trận 𝐶 𝑘 (1.16) và các véc tơ v ế phải 𝐹𝑗 𝑘 được tính theo các công thức truy toán 𝐶 𝑘 = 𝐶 𝑘−1 2 𝑘−1 − 2𝐸; 𝐹𝑗 𝑘 = 𝐹𝑗 −2 𝑘−1 + 𝐶 𝑘−1 𝑘−1 𝐹𝑗 𝑘−1 + 𝐹𝑗 +2 𝑘−1 , 𝑗 = 2𝑘 , 2. 2𝑘 , 3. 2𝑘 , … , 𝑁 − 2𝑘 , 𝑘 = 1, 2, 3, … (1.17) Từ (1.15) suy ra rằng sau 𝑛 − 1 bước khử (𝑙 = 𝑛 − 1) ta thu được hệ chỉ gồm một phương trình đối với biến 𝑌2𝑁 −1 = 𝑌𝑁/2 là: 𝐶 𝑛 −1 𝑌𝑗 = 𝐹𝑗 𝑛−1 + 𝑌𝑗 −2𝑛 −1 + 𝑌𝑗 +2𝑛 −1 = 𝐹𝑗 𝑛−1 + 𝑌0 + 𝑌𝑁 ; 𝑌0 = 𝐹0 , 𝑌𝑁 = 𝐹𝑁 (1.18) với vế phải đã biết. Vì vậy từ (1.18) ta có thể tìm được 𝑌𝑁/2 , kết hợp với (1.15) ta thấy rằng tất cả các ẩn được tìm liên tiếp từ các phương trình 𝐶 𝑘−1 𝑌𝑗 = 𝐹𝑗 𝑘−1 + 𝑌𝑗 −2𝑘−1 + 𝑌𝑗 +2𝑘−1 , 𝑌0 = 𝐹0 , 𝑌𝑁 = 𝐹𝑁 , 𝑗 = 2𝑘−1 , 3. 2𝑘−1 , 5. 2𝑘−1 , … , 𝑁 − 2𝑘−1 , 𝑘 = 𝑛, 𝑛 − 1, … , 1 Việc tính các 𝐹𝑗 𝑘 (1.19) theo công thức truy toán (1.17) có thể dẫn đến việc tích lũy sai số nếu như chuẩn của ma trận 𝐶 𝑘−1 lớn hơn 1. Ngoài ra các ma trận 𝐶 𝑘 nói chung là các ma trận đầy đủ. Điều này dẫn đến việc tăng khối lượng tính toán khi tính các 𝐹𝑗 𝑘 theo (1.17). Để khắc phục những khó khăn trên, thay cho 𝐹𝑗 𝑘 ta sẽ tính các véc tơ 𝑝𝑗𝑘 và 𝑞𝑗𝑘 liên hệ với 𝐹𝑗 𝑘 theo công thức sau: Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn 14 𝐹𝑗 𝑘 𝑘 =𝐶 𝑝𝑗𝑘 + 𝑞𝑗𝑘 , 𝑗 = 2𝑘 , 2. 2𝑘 , 3. 2𝑘 , … , 𝑁 − 2𝑘 , 𝑘 = 0, 1, 2, 3, … , 𝑛 − 1 (1.20) (0) trong đó ta chọn 𝑝𝑗 (0) = 0 và 𝑝𝑗 = 𝐹𝑗 , 𝑗 = 1, 2, … , 𝑁 − 1. Để tìm mối quan hệ giữa 𝑝𝑗𝑘 và 𝑞𝑗𝑘 ta thay (1.20) vào (1.17) ta được: 𝐶 (𝑘) 𝑘 (𝑘) 𝑝𝑗 + 𝑞𝑗 =𝐶 𝑘+1 𝑞𝑗 𝑘−1 𝑘−1 + 𝑝𝑗 −2𝑘−1 + 𝐶 𝑘−1 𝑝𝑗 𝑘−1 𝑘−1 + 𝑝𝑗 +2𝑘−1 + 𝑘−1 𝑘−1 𝑘 𝑘 𝑘 𝑘 𝑞𝑗 −2 𝑘−1 + 𝑞𝑗 +2𝑘−1 , 𝑗 = 2 , 2. 2 , 3. 2 , … , 𝑁 − 2 , 𝑘 = 0, 1, 2, … , 𝑛 − 1 (1.21) (𝑘) (𝑘) Ta sẽ chọn 𝑝𝑗 (𝑘) 𝑞𝑗 và 𝑞𝑗 (𝑘) = 2𝑝𝑗 thỏa mãn 𝑘−1 𝑘−1 + 𝑞𝑗 −2 𝑘−1 + 𝑞𝑗 +2𝑘−1 Khi đó, kết hợp với công thức 𝐶 𝐶 (𝑘) 𝑘 (𝑘−1) = 𝑝𝑗 𝑘−1 (𝑘−1) + 𝑝𝑗 (𝑘−1) 𝑆𝑗 + 2𝐸 = 𝐶 𝑘−1 + 𝑝𝑗 −2 𝑘−1 + 𝐶 = 𝑞𝑗 (𝑘) Đặt 𝑆𝑗 𝐶 (𝑘−1) 𝑝𝑗 𝑘 (1.22) 𝑘−1 (𝑘−1) (𝑘−1) ta có: 𝑘−1 𝑝𝑗 𝑘−1 + 𝑝𝑗 +2 𝑘−1 từ (1.23) suy ra 𝑆𝑗 = 𝑞𝑗 2 𝑘+1 (1.23) phải thỏa mãn 𝑘−1 𝑘−1 + 𝑝𝑗 −2 𝑘−1 + 𝑝𝑗 +2𝑘−1 (1.24) Kết hợp (1.20) và (1.24) ta có thuật toán xác định các véc tơ 𝑝𝑗𝑘 và 𝑞𝑗𝑘 𝐶 𝑘−1 (𝑘) 𝑝𝑗 (𝑘−1) 𝑆𝑗 (𝑘−1) = 𝑞𝑗 (𝑘−1) = 𝑝𝑗 (𝑘−1) + 𝑆𝑗 (0) 𝑘−1 𝑘−1 + 𝑝𝑗 −2 𝑘−1 + 𝑝𝑗 +2𝑘−1 , 𝑞𝑗 (𝑘) ; 𝑞𝑗 (𝑘) = 2𝑝𝑗 (0) = 𝐹𝑗 , 𝑝𝑗 = 0, 𝑘−1 𝑘−1 + 𝑞𝑗 −2 𝑘−1 + 𝑞𝑗 +2𝑘−1 ; (1.25) 𝑗 = 2𝑘 , 2. 2𝑘 , 3. 2𝑘 , … , 𝑁 − 2𝑘 , 𝑘 = 0, 1, 2, … , 𝑛 − 1 (𝑘−1) Đặt 𝑡𝑗 (𝑘−1) = 𝑌𝑗 − 𝑝𝑗 , thay công thức (1.20) vào (1.19) ta thấy 𝑌𝑗 có thể tính được từ các công thức sau: 𝐶 𝑘−1 (𝑘−1) 𝑡𝑗 (𝑘−1) = 𝑞𝑗 (𝑘−1) + 𝑌𝑗 −2𝑘−1 + 𝑌𝑗 +2𝑘−1 , 𝑌𝑗 = 𝑝𝑗 + 𝑡𝑗𝑘−1 𝑌0 = 𝐹0 , 𝑌𝑁 = 𝐹𝑁 , 𝑗 = 2𝑘 , 2. 2𝑘 , 3. 2𝑘 , … , 𝑁 − 2𝑘 , 𝑘 = 𝑛, 𝑛 − 1, … , 1 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên (1.26) http://www.lrc-tnu.edu.vn 15