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 ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC CNTT & TT
BÙI THỊ THU TRANG
MỘT SỐ PHƯƠNG PHÁP XÁC ĐỊNH NGHIỆM GẦN ĐÚNG
ĐỐI VỚI BÀI TOÁN BIÊN ELLIPTIC CẤP HAI TRONG
MIỀN PHỨC TẠP HOẶC ĐIỀU KIỆN BIÊN PHỨC TẠP LUẬN VĂN THẠC SỸ KHOA HỌC MÁY TÍNH Chuyên ngành : Khoa học máy tính
Mã số : 60 48 01
1.4.2. Bài toán biên Neumann 19
1.5. Phƣơng pháp lặp giải phƣơng trình toán tử 24
1.5.1. Lƣợc đồ lặp hai lớp 24
1.5.2. Lƣợc đồ dừng, định lý cơ bản về sự hội tụ của phép lặp . 27
Chƣơng 2 30
CÁC PHƢƠNG PHÁP XẤP XỈ BIÊN 30
ĐỐI VỚI BÀI TOÁN CÓ ĐIỀU KIỆN BIÊN KÌ DỊ 30
2.1. Cơ sở của phƣơng pháp 30
2.2. Các phƣơng pháp xấp xỉ biên (BAMs) 31
2.2.1. Cơ sở phƣơng pháp 31
2.2.2. Các phƣơng pháp BAMs 32
2.3. Phƣơng pháp xấp xỉ (GFIFs) 32
2.4. Giới thiệu bài toán Motz: 34
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
ii
2.4.1. Kết quả sử dụng các phƣơng pháp BAMs 34
2.4.2. Kết quả sử dụng phƣơng pháp GFIFs 37
Chƣơng 3 41
PHƢƠNG PHÁP CHIA MIỀN 41
GIẢI BÀI TOÁN BIÊN VỚI BIÊN KỲ DỊ 41
3.1. Cơ sở của phƣơng pháp 41
3.2. Ứng dụng của phƣơng pháp chia miền đối với bài toán Motz 43
3.3. Mở rộng phƣơng pháp chia miền trong trƣờng hợp tổng quát: 47
PHẦN KẾT LUẬN 52
TÀI LIỆU THAM KHẢO 53
PHẦN PHỤ LỤC 55
k
39
Hình 3.1. Mô tả phƣơng pháp chia miền 41
Hình 3.2. Mô hình Bài toán Mozt áp dụng phƣơng pháp chia miền 44
Hình 3.3. Nghiệm của bài toán Motz với phƣơng pháp chia miền. 46
Hình 3.4. Mô hình bài toán Motz dạng tổng quát 47
Hình 3.5. Đồ thị nghiệm trong trƣờng hợp tổng quát. 50
Hình 3.6. Đƣờng cong đạo hàm bậc nhất 51
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
MỞ ĐẦU
Một số mô hình trong vật lý và cơ học đều đƣa đến việc xác định nghiệm của
bài toán biên elliptic cấp hai với các điều kiện biên khác nhau. Trong những trƣờng
hợp đơn giản khi miền hình học là miền hình chữ nhật và có điều kiện biên hỗn hợp
yếu (tức là trên một cạnh chỉ có 1 loại điều kiện biên Dirichlet hoặc Neumann) thì
bằng việc sử dụng các phƣơng pháp sai phân, bài toán vi phân đƣợc đƣa về hệ
phƣơng trình véc tơ ba điểm và xác định nghiệm gần đúng qua việc giải các hệ
phƣơng trình đại số tuyền tính bằng các thuật toán đã biết. Tuy nhiên trong miền
hình học không phải là miền hình chữ nhật hoặc điều kiện biên là hỗn hợp mạnh
(tức là trên một cạnh có hai loại điều kiện biên Dirichlet và Neumann) thì bài toán
xuất hiện điểm kì dị tại các góc của miền hoặc điểm phân chia giữa hai loại điều
kiện biên khi đó phƣơng pháp sai phân sẽ gặp khó khăn. Để giải quyết các dạng bài
toán trên, chúng ta có thể sử dụng phƣơng pháp chia miền tìm nghiệm xấp xỉ bằng
tổng hữu hạn các hàm cơ sở xung quanh điểm kì dị. Nhƣ vậy, việc tìm hiểu, nghiên
cứu các phƣơng pháp tìm nghiệm xấp xỉ của các bài toán biên elliptic cấp hai trong
trƣờng hợp miền phức tạp hoặc điều kiện biên phức tạp cũng nhƣ so sánh giữa các
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
3
Chƣơng 1
MỘT SỐ KIẾN THỨC CƠ BẢN VỀ GIẢI SỐ
PHƢƠNG TRÌNH ĐẠO HÀM RIÊNG
Trong chƣơng này, luận văn trình bày một số kiến thức liên quan đến việc
giải số phƣơng trình đạo hàm riêng bao gồm: Cơ sở về phƣơng pháp lƣới, phƣơng
pháp lặp, thuật toán thu gọn khối lƣợng tính toán, giới thiệu các thƣ viện chƣơng
trình TK2004 và RC2009 tìm nghiệm số của các bài toán biên hỗn hợp trong miền
chữ nhật. Những kiến thức cơ sở và kết quả đƣợc tham khảo từ các tài liệu
[3,4,8,10].
1.1. Phƣơng pháp sai phân.
1.1.1. Lƣới sai phân.
Xét bài toán
u f , x ,
u g, x .
(1.1)
trong đó
( x,y ) R ,a x b,c y d
, tập
hkhkhk
=
gọi là một
lƣới sai phân trên
.
1.1.2. Hàm lƣới.
Mỗi hàm số xác định tại các nút của lƣới gọi là một hàm lƣới, giá trị của hàm
lƣới
),( yxu
tại nút lƣới
( , )ij
viết tắt là
i , j
u
,
MjNi ,0,,0
. Mỗi hàm
( , )u i j
xác định tại mọi
),( yx
tạo ra hàm lƣới
u
xác định bởi
i , j
u
.
1
4
4
),(
,
constCyx
y
u
max
yx
=),(
2
4
),(
.
Do đó theo công thức Taylor ta có:
2 2 3 3
4
1
23
( , ) = ( , ) = ( , ) ( )
Một cách tƣơng tự:
i j i j
u( x ,y ) u( x ,y k )
1
=
2 2 3 3
4
23
= ( , ) ( )
2! 3!
ij
u k u k u
u x y k k
y y y
i j i j
u( x ,y ) u( x ,y k )
1
2 2 3 3
4
Vậy ta có:
11
2
( , ) 2 ( , ) ( , )
i j i j i j
u x y u x y u x y
h
11
22
2
( , ) 2 ( , ) ( , )
()
i j i j i j
u x y u x y u x y
u h k
k
Ta đặt:
1 , 1, , 1 , , 1
kh
xấp xỉ
toán tử
, điều đó cho phép thay phƣơng trình vi phân bằng phƣơng trình sai phân:
hk ij ij i j i j hk
u f , f f ( x ,y ), ( x ,y ) ==
tức là:
1, , 1, , 1 , , 1
22
22
( , ), ( , )
i j i j i j i j i j i j
i j i j hk
u u u u u u
f x y x y
hk
(1.2)
đồng thời thay điều kiện biên bằng điều kiện:
ij i j i j hk
u g( x ,y ), ( x ,y )
(1.3)
Y
là véc tơ cần tìm,
C
là ma trận vuông,
j
F
là véc tơ cho trƣớc. Ý
tƣởng của phƣơng pháp rút gọn hoàn toàn giải (1.1) là khử liên tiếp các ẩn
j
Y
đầu
tiên với các
j
lẻ, sau đó từ các phƣơng trình còn lại khử các
j
Y
với
j
là bội của 2,
rồi 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
n
N 2=
thì sau một số lần khử sẽ còn lại một phƣơng trình chứa véc tơ ẩn
/2N
Y
mà từ đó
/2N
Y
6
rút gọn hoàn toàn là một biến thể của phƣơng pháp khử Gauss áp dụng cho bài toán
(1.4) trong đó việc khử các biến đƣợc thực hiện theo một thứ tự đặc biệt. Sau đây, ta
sẽ mô tả cụ thể phƣơng pháp.
Giả sử
0>,2= nN
n
. Ký hiệu
11,2, ,=;=,=
(0)(0)
NjFFCC
jj
. Khi đó (1.4)
đƣợc viết dƣới dạng
1)(1=
(0)
1
0
1
NjFYYCY
jjjj
,
00
= FY
,
NN
FY =
(0)
=
jjjj
FYYCY
Nhân 2 vế của phƣơng trình thứ hai với
(0)
C
vào bên trái rồi cộng cả 3 phƣơng trình
lại ta đƣợc:
22,4, ,=,=
(1)
12
(1)
2
NjFYYCY
jjjj
,
00
= FY
,
NN
FY =
(1.6)
trong đó:
j
Y
với
j
lẻ sẽ tìm đƣợc từ phƣơng trình
11,3, ,=,=
11
(0)(0)
NjYYFYC
jjjj
(1.7)
Nhƣ vậy hệ (1.5) tƣơng đƣơng với hệ gồm (1.6) và (1.7).
Bƣớc khử thứ hai: ở bƣớc khử này ta sẽ tiến hành khử các
j
Y
của hệ (1.6)
với
j
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.6)
(1)
22
(1)
4
=
trình lại ta đƣợc
44,8, ,=,=
(2)
4
(2)
4
NjFYYCY
jjjj
,
00
= FY
,
NN
FY =
(1.8)
trong đó:
ECC 2)(=
2)(1(2)
,
44,8, ,=,=
(1)
2
(1)(1)(1)
2
(2)
.
Cứ tiếp tục quá trình khử này. Kết quả là sau bƣớc khử thứ
l
ta nhận đƣợc
một hệ gồm
1
l
C
N
ẩn trong đó
j
là bội của
l
2lllll
jl
j
j
l
l
j
NjFYYCY 2, ,,3.2,2.22=,=
)(
2
)(
2
111
Nj
kkk
,
1, ,1,,= llk
(1.10)
trong đó các ma trận
)(k
C
và các véc tơ vế phải
)(k
j
F
đƣợc tính theo các công thức
truy toán:
ECC
kk
2)(=
21)()(
,
1
2
1)(11)(
1
2
)(
=1
N
N
YY
là
NNN
n
jn
j
n
j
n
jj
n
FYFYYYFYYFYC ==,==
000
1)(
1
2
1
2
1)(1)(
(1.12)
(1.13)
Các công thức trên đã mô tả phƣơng pháp rút gọn hoàn toàn. Việc tính các
)(k
j
F
theo công thức truy toán có thể dẫn đến việc tích luỹ sai số nếu nhƣ chuẩn của
ma trận
1)( k
C
lớn hơn 1. Ngoài ra các ma trận
)(k
C
nói chung là các ma trận đầy đủ,
thậm chí cả với ma trận ban đầu là
CC =
(0)
(ma trận ba đƣờng chéo). Điều này dẫn
đến tăng khối lƣợng tính toán khi tính các
)(k
j
F
theo (1.13). Để khắc phục những
khó khăn trên, thay cho
jj
. Bằng các công thức toán học, có
thể thấy mối quan hệ mà
)(k
j
p
và
)(k
j
q
phải thoả mãn nhƣ sau.
=
)()()( k
j
k
j
k
qpC
1)(
1
2
1)(
1
2
1)(
1
2
1)(1)(1)(
1
j
k
j
k
qqppCpqC
,
1,2,3 =,2, ,,3.2,2.22= kNj
kkkk
,
Ta sẽ chọn
)(k
j
p
và
)(k
j
q
thoả mãn
1)(
1)(2
1)(
1)(2
)()(
2=
k
k
j
k
j
kk
k
j
k
j
k
j
k
ppCpqpC
Đặt
1)()(1)(
=
k
j
k
j
k
j
k
j
k
ppqSC
Nhƣ vậy ta thu đƣợc thuật toán sau đây để xác định các véc tơ
)(k
j
p
và
)(k
j
q
1)(
1
2
1)(
1
2
1)(1)(1)(
=
2=
k
kj
k
kj
k
j
k
j
qqpq
(1.15)
0=;=
(0)(0)
jjj
pFq
,
kkkk
Nj 2, ,,3.2,2.22=
,
10,1,2, ,= nk
.
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
k
j
k
j
k
YYqtC
,
1)(1)(
=
k
j
k
jj
tpY
,
NN
FYFY =;=
00
,
kkkk
Nj 2, ,,3.2,2.22=
,
1, ,1,= nnk
(1.16)
Nhận xét rằng các quá trình (1.15) và (1.16) luôn cần tính ma trận nghịch đảo
1)1(
][
k
CC
, trong đó
E
l
CC
k
kl
2
1)(2
cos2=
1,
.
Nhƣ vậy, chẳng hạn ta có phƣơng trình
=
1)( k
C
, (1.17)
thì với việc giải lần lƣợt các phƣơng trình
1
11,
1,2, ,2=,=
jj
qCp
và tính
22,4,6, ,=(0),
1)(
(0)
1)(
(1)(1)
2=
Nj
jjjj
qqpq
Bƣớc 3: Với
12,3, ,= nk
xác định các véc tơ
kkkkk
k
j
k
k
j
k
jj
Njppq 2 ,,3.2,2.22=,=
1)(
1
l
j
l
jkl
C
Khi đó
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
(21)()(
=
k
j
k
j
k
j
pp
,
kkk
k
j
k
j
k
jj
NjYYq
Sau đó, với mỗi
1
1,2, ,2=
k
l
và với mỗi
111
2 ,3.22=
kkk
Nj
, giải phƣơng
trình
11,3,5, ,=,=
11
(0)
NjYYqCY
jjjj
1.2.2. Bài toán biên thứ hai
Xét bài toán thứ hai
.2
,11,
,0,
1
11
00
NNN
jjjj
FCYY
nnkNjYYFYC
kkk
jj
k
jj
k
kk
Trong đó
)(k
j
F
và
)(k
C
đƣợc xác định bởi công thức truy toán 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
11
,2][
,2
2, ,2.2,2
,
,
2)1()(
)1()1(1
j
k
j
k
k
kk
(1.20)
Kí hiệu:
,
)()()()( k
j
k
j
kk
j
)1()1()1(
)1()1(
k
j
k
j
k
j
k
j
k
kk
qpqSC
,
)1()1()(
k
j
k
j
k
j
jjj
kkk
pFq
nkNj
Tƣơng tự, với j=N, ta có:
,0;
,22
,
,2
)0()0(
)1(
2
)()(
)1()1()(
)1(
2
)1(1)1(
)1(
)1(
Trong đó
1,2, ,1,,2, ,2.3,2,
111)1()1(
nnkNjtpY
kkkk
j
k
jj
Trong đó
)1( k
j
t
là nghiệm của phƣơng trình.
)1()1(
22
)1()1()1(
kk
jj
k
j
k
j
)1()0(
)1()1(
Sau đó, với
)1(
2, ,2,1
k
l
và với mỗi
,2, ,2.2,2
kkk
Nj
giải phƣơng trình
)1()(
1,
l
j
l
Bƣớc 3: Với k=1,2,…,n-1 xác định các véc tơ
)1(
2
)1()0(
)1(
2
k
N
k
NN
k
pqv
Sau đó, với
,2, ,2,1
)1(
k
l
giải phƣơng trình
k
N
k
N
k
k
qpqvpP
Quá trình ngƣợc:
Bƣớc 1: Xác định Y
N
, xác định véc tơ
0
)()0(
2Yqv
n
NN
sau đó, với l=1,2,…,n, giải hệ
)1()(
,
l
N
l
Nnl
vvC
NjYYqv
kk
Sau đó, với
1
2, ,2,1
k
l
và với mỗi
,2, ,2.2,2
111
kkk
Nj
giải phƣơng
trình
)1()(
1,
l
j
l
jkl
vvC
Khi đó
,.
u cu f x
ux
(1.21)
trong đó
22
( ); ( )f L L
Rời rạc hoá miền
bằng lƣới
h
1 2 1 1 2 2
= = ( , ) = ( , ); = ; = ;
h i j i j
x i j x x x ih x jhSố hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên http://www.lrc-tnu.edu.vn
= ( ; ; ; )
= ( ; ; ; )
= ( ; ; ; )
j j j j
NN
T
j j j M j
T
M
T
N N N M N
Y CY Y F j N
Y F Y F
Y y y y
F g g g
F g g g
2
2 1, 0,
2
2 2,
2
2 2,
2
2 1, ,
jj
j
j
Mj
M j M j
h rg
h
F
h
h rg
0000( , 1, 2, 3, 4, 1, 2, , , , , 1, 2, 1, 2)u phi b b b b L L c M N n p p q q
trả lại
nghiệm bằng số của bài toán, trong đó:
+
MN
là lƣới trên hình chữ nhật.
+
phi
là ma trận giá trị của hàm vế phải.
+
1, 2, 3, 4b b b b
là véc tơ giá trị của điều kiện biên trên các cạnh.
+
1, 2LL
là kích thƣớc của hình chữ nhật.
+
1, 2, 1, 2p p q q
là tọa độ của các đỉnh hình chữ nhật trong không gian lƣới.
1.3.2. Bài toán biên hỗn hợp
Cho miền
1 2 1 1 2 2
= = ( , );0 < < ;0 < <x x x x l x l
Xét bài toán biên hỗn hợp
1
2
2
= ( );
= ( );
1 1 1 1 1 2 2 2
= = ( ,0);0 = ( , );0x x x l x l x x l
2 2 2
= (0, );0x x x l
2 1 2 1 1
= = ( , );0x x l x l
Tƣơng tự, ta đƣa bài toán sai phân về hệ phƣơng trình vec tơ 3 điểm
00
11
1
=
= 1 1
2 = =
j j j j
N N N
YF
Y CY Y F j N
Y CY F j N
rr
r r r
rr
rr
r r r
rr
2
2 1, 0,
2
2 2,
2
2 1, 2 1, 0,
2
2 2, 2 2,
2
2 2, 2 2,
2
2 1, 2 1, ,
2
2
2
2
N N N
NN
N
M N M N
M N M N M N
h h rg
hh
F
hh
h h rg
là véc tơ giá trị của điều kiện biên trên các cạnh.
+
1, 2LL
là kích thƣớc của hình chữ nhật.
+
1, 2, 1, 2p p q q
là tọa độ các đỉnh của hình chữ nhật trong không gian 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
17
Hoàn toàn tƣơng tự, trong các trƣờng hợp khi điều kiện biên Neumann trên 2
cạnh, 3 cạnh hoặc 4 cạnh của hình chữ nhật, sử dụng phƣơng pháp sai phân, chúng
ta đƣa bài toán elliptic về các hệ phƣơng trình vecto 3 điểm và sử dụng thuật toán
thu gọn khối lƣợng tính toán xác định nghiệm bằng số trên từng điểm lƣới. Các hàm
đƣợc thiết lập tƣơng tự là
0010( ), 0100( ), 1000( ), 1100( ), , 1111( )u u u u u
sẽ cho
phép xác định nghiệm số với các bài toán tƣơng ứng. Trong các kí hiệu hàm trên, kí
hiệu 0 tƣơng ứng với biên Dirichlet, kí hiệu 1 tƣơng ứng với biên Neumann. Các hệ
thống hàm trên đã đƣợc xây dựng thành thƣ viện TK2004 cho phép giải số tất cả
các bài toán biên elliptic với điều kiện biên hỗn hợp, các kết quả trên đã đƣợc đƣa ra
trong [3].
1.4. Giới thiệu thƣ viện RC2009
Thƣ viện chƣơng trình RC2009 là sự phát triển của thƣ viên TK2004 tìm
nghiệm số của bài toán biên hỗn hợp trong trƣờng hợp toán tử của phƣơng trình
phức tạp hơn.
Cho
là hình chữ nhật
.
1.4.1. Bài toán biên Dirichlet
Xét trƣờng hợp khi toán tử
lu u
tức là điều kiện biên dạng Dirichlet,
12
k ,k ,c
là các hằng số,
là hình chữ nhật có kích thƣớc hai cạnh là L
1
, L
2
. Xuất phát từ
phƣơng pháp lƣới, chia miền
thành
MN
điểm lƣới, trong đó
2 , 0
n
Nn
. Kí
hiệu
12
12
là các vé c tơ cấp (M-1), C là ma trận hệ số cấp
11MM
đƣợc xác định nhƣ sau
1, 2, 1,
0 1,0 2,0 1,0 1, 2, 1,
, , , , 0,
, , , , , , ,
j j j M j
M N N N M N
Y u u u j N
F g g g F g g g
Ma trận C có dạng
0 0 0 0
0 0 0
0 0 0 0
,
0 0 0 0
0 0 0
0 0 0 0
dr
r d r
rd
2
2,
2
2
2
2,
2
2
2
1, ,
2
jj
j
j
Mj
M j M j
h
rg
k
h
k
F
h
k
h
rg
k
, 2 1
k h h
r d r c
k h k
.
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
19
Trên cơ sở thuật toán thứ nhất tiến hành cài đặt giải hệ phƣơng trình trên.
Thiết kế các hàm RC0000(phi,b1,b2,b3,b4,l1,l2,k1,k2,C,N,M,n), thực hiện thuật
toán thu gọn.
Hàm v0000(phi,b1,b2,b3,b4,l1,l2,k1,k2,C,M,N,n,p1,p2,q1,q2) trả lại ma trận
nghiệm xấp xỉ của bài toán (1.24) bắt đầu từ tọa độ (p
1
,q
1
) đến (p
2
,q
2
).
1.4.2. Bài toán biên Neumann
Xét bài toán biên hỗn hợp
22
12
22
12
Neumann.
Từ phƣơng pháp sai phân với độ chính xác
22
12
O h h
chuyển bài toán vi
phân (1.25) về bài toán sai phân tƣơng ứng với hệ phƣơng trình vé c tơ ba điểm
1 1 0 0 1
; , 2 ; 1, 1
j j j j N N N
Y CY Y F Y F Y CY F j N
trong đó Y
j
là các véc tơ nghiệm,
j
F
là các vé c tơ cấp
1M
,C là ma trận hệ số cấp
11MM
đƣợc xác định nhƣ sau:
2
2
2
1, 2 4 2
2
21
22
22
21
N
N
N
MN
MN
h
h b rb N
k
h
hb
k
F
h
h b M
k
h
h b M rb N
k
0 0 0
0 0 0 0
dr
r d r
rd
C
dr
r d r
rd
h
rb j
k
1,2, 1jN