Luận văn Thạc sĩ Toán học: Về véctơ trọng số RBF cho phương pháp không lưới RBF – FD trong không gian ba chiều

ĐẠI HỌC THÁI NGUYÊN TRƯỜNG ĐẠI HỌC KHOA HỌC ————— o0o —————

VŨ XUÂN HIỂN

VỀ VÉCTƠ TRỌNG SỐ RBF CHO

PHƯƠNG PHÁP KHÔNG LƯỚI RBF-FD TRONG KHÔNG GIAN BA CHIỀU

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

THÁI NGUYÊN, 6/2020

ĐẠI HỌC THÁI NGUYÊN TRƯỜNG ĐẠI HỌC KHOA HỌC ————— o0o —————

VŨ XUÂN HIỂN

VỀ VÉCTƠ TRỌNG SỐ RBF CHO

PHƯƠNG PHÁP KHÔNG LƯỚI RBF-FD TRONG KHÔNG GIAN BA CHIỀU

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

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

NGƯỜI HƯỚNG DẪN KHOA HỌC:

TS. ĐẶNG THỊ OANH

Thái Nguyên, 6/2020

ii

Mục lục

iv Lời cảm ơn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Danh mục ký hiệu . . . . . . . . . . . . . . . . . . . . . . . . v

Mở đầu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

Chương 1. Một số kiến thức chuẩn bị 3

1.1 Cơ sở của bài toán nội suy hàm số với dữ liệu phân tán . . 3 1.2 Một số định nghĩa và khái niệm . . . . . . . . . . . . . . . 4

1.3 Lược đồ phương pháp sai phân hữu hạn giải phương trình

đạo hàm riêng trong không gian 3 chiều . . . . . . . . . . 10

1.4 Kết luận . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Chương 2. Phương pháp không lưới RBF – FD giải phương

bố không đều . . . . . . . . . . . . . . . . . . . . 16 2.2.2 Nội suy không có thành phần đa thức . . . . . . 18 2.2.3 Nội suy có thành phần đa thức . . . . . . . . . . 20 2.2.4 Véctơ trọng số RBF trong không gian ba chiều . . . 23 . . . . . . . . 26

15 trình Poisson trong không gian ba chiều 2.1 Bài toán mở đầu . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Véctơ trọng số từ nội suy hàm cơ sở bán kính . . . . . 16 2.2.1 Véctơ trọng số từ vi phân số trên các tâm phân

2.3 Thuật toán chọn tâm dựa trên các góc khối

2.4

Lược đồ RBF-FD giải bài toán Elliptic trong không gian ba chiều . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.5 Thử nghiệm số . . . . . . . . . . . . . . . . . . . . . . . . 27

Kết luận 31

iii

Tài liệu tham khảo 32

iv

Lời cảm ơn

Luận văn được hoàn thành tại Trường Đại học Khoa học - Đại học

Thái Nguyên. Trong quá trình học tập và thực hiện luận văn này, Trường Đại học Khoa học đã tạo mọi điều kiện tốt nhất để tác giả học tập,

nghiên cứu. Tác giả xin được bày tỏ lòng biết ơn chân thành đến các thầy, cô trong khoa Toán - Tin, trong Trường Đại học Khoa học - Đại

học Thái Nguyên. Đặc biệt, tác giả xin bày tỏ lòng biết ơn sâu sắc tới TS. Đặng Thị Oanh - Người đã tận tình hướng dẫn tác giả hoàn thành

luận văn này. Tác giả cũng xin được gửi lời cảm ơn tới Ban giám hiệu trường THPT Ân Thi, Hưng Yên và tập thể các thầy cô giáo trong tổ

Toán Tin của Trường đã tạo điều kiện giúp đỡ tác giả trong thời gian tác giả tham gia học cao học.

Cuối cùng tôi xin gửi lời cám ơn đến tập thể lớp K12A6, gia đình bạn

bè đã giúp đỡ, động viên tôi trong suốt quá trình học tập và thực hiện luận văn này.

Thái Nguyên, tháng 6 năm 2020 Tác giả

Vũ Xuân Hiển

v

Danh mục ký hiệu

Miền hình học. Ω

Tập các các tâm trong miền và trên biên Ω. Tập các tâm nằm trong miền Ω.

Ξ Ξint ∂Ω ζ

g f Tập các tâm nằm trên biên ∂Ω. Tâm thuộc tập Ξint. Hàm trên biên. Hàm vế phải của phương trình Poisson.

w u véc tơ trọng số. Nghiệm giải tích.

˜u Nghiệm xấp xỉ.

∞ Rn Vô cùng. Không gian n chiều.

λ φ Giá trị riêng của ma trận. Hàm cơ sở bán kính.

(cid:15) A Tham số hình dạng. Ma trận của hệ phương trình đại số tuyến tính.

b x Véc tơ vế phải của hệ phương trình đại số tuyến tính. Nghiệm của hệ phương trình đại số tuyến tính.

E X

k Ma trận đơn vị. Bộ tâm gồm ξ và ζ. Ký hiệu: X = {ζ, ξ1, ..., ξk} . Số các điểm ξi cần thiết trong tập Ξζ.

1

Mở đầu

Trong suốt thế kỷ XX, một loạt các phương pháp số đã hình thành

và phát triển như các phương pháp sai phân hữu hạn (finite difference - FD), phương pháp phần tử hữu hạn (finite element method - FEM),

phương pháp thể tích hữu hạn (Finite Volume Method -FVM), phương pháp phần tử biên (Boundary Element Method - BEM). . . Những phương

pháp này đã đem lại những đóng góp to lớn trong việc ứng dụng toán học vào thực tiễn. Tuy nhiên, chúng còn nhiều hạn chế khi áp dụng vào

lớp các bài toán thực tế có cấu trúc phức tạp như: lưới biến dạng trên phạm vi rộng, số chiều không gian cao, hàm vế phải hoặc hàm điều kiện

biên có kì dị (có độ dao động lớn). Khó khăn lớn nhất ở đây là sinh lưới, duy trì lưới và cập nhật lưới. Đó là lý do thúc đẩy các nhà khoa học

thuộc các lĩnh vực khác nhau, tìm kiếm những phương pháp mới nhằm

khắc phục những hạn chế này của các phương pháp lưới.

Để khắc phục một số nhược điểm của phương pháp lưới, người ta

đã đưa ra phương pháp không lưới giải phương trình đạo hàm riêng [4, 5, 9, 11, 12, 14]. Một trong các cách tiếp cận không lưới là phương

pháp RBF – FD (Radial Basis Funcion - Finite Difference) [9, 11, 12, 14]. - Hàm cơ sở bán kính Φ : Rd → R là hàm xác định dương và giá trị của nó chỉ phụ thuộc vào chuẩn của véctơ biến, tức là Φ(x) = φ(||x||2) với mọi x ∈ Rd và φ : [0; ∞) → R là một hàm cho trước nào đó [6, 10]. Sử dụng sự thay đổi của hàm này để xây dựng nội suy hàm cơ sở bán kính (RBF – Interpolation) [6, 10]. Tiếp đó nội suy hàm cơ sở bán kính được

sử dụng xấp xỉ toán tử vi phân, nhằm tạo ra các phương pháp xấp xỉ

giải phương trình vi phân đạo hàm riêng. Phương pháp RBF – FD được xây dựng theo lược đồ này. Phương pháp RBF – FD là phương pháp

không lưới sử dụng nội suy hàm cơ sở bán kính với cách tiếp cận địa phương, và dựa trên sự rời rạc hóa giống như phương pháp FD để tính

xấp xỉ nghiệm tại một số điểm rời rạc trong miền xác định. Khi sử dụng phương pháp RBF – FD giải bài toán d chiều với d lớn tùy ý, thì thay

vì phải làm việc với hàm d biến, ta chỉ cần làm việc với hàm một biến. Một lợi thế của kỹ thuật rời rạc không lưới là chỉ cần dựa trên tập điểm

2

độc lập phân bố bất kỳ, không cần tạo ra cấu trúc lưới. Do đó, không còn cần chi phí dành cho sinh lưới, duy trì lưới và cập nhật lưới.

- Các nghiên cứu về phương pháp không lưới dựa trên hàm cơ sở bán kính (Radial Basis Function – RBF) và ứng dụng giải phương trình đạo

hàm riêng đã được nhiều nhà khoa học trong và ngoài nước quan tâm. Luận văn tập trung trình bày cách tính véctơ trọng số dựa trên nội

suy RBF trong không gian ba chiều cho phương pháp sai phân hữu hạn

không lưới trên cơ sở đọc hiểu và tổng hợp từ các quyển sách [1, 6, 7, 10] và các bài báo [14, 15].

3

Chương 1

Một số kiến thức chuẩn bị

Trong chương này chúng tôi trình bày một số kiến thức cơ bản về cơ sở của bài toán nội suy hàm số với dữ liệu phân tán; nội suy hàm số với

dữ liệu phân tán; một số hàm cơ sở bán kính và lược đồ phương pháp

sai phân hữu hạn giải phương trình đạo hàm riêng trong không gian 3 chiều. Nội dung chính của chương này được tham khảo chủ yếu từ các

1.1 Cơ sở của bài toán nội suy hàm số với dữ liệu

phân tán

tài liệu [1, 6, 7, 10].

Trong thực tế nhiều khi phải phục hồi một hàm số f (x) tại mọi giá

trị của x trên đoạn [a; b] mà chỉ biết một số hữu hạn giá trị của hàm số tại một số hữu hạn các điểm rời rạc của đoạn đó. Các giá trị đó được

cung cấp qua thực nghiệm hay tính toán, vì vậy nảy sinh một vấn đề toán học như sau: Trên đoạn [a; b] cho một lưới các điểm chia (điểm nút) xi, i = 0, 1, 2, · · · , n và tại các nút xi cho giá trị của hàm số y = f (x) là yi = f (xi), i = 0, 1, 2, · · · , n. Cần xây dựng đa thức nội suy Pn(x) sao cho Pn(x) trùng với f (x) tại tất cả các nút xi, nghĩa là:

i = 0, 1, 2, · · · , n. Pn(xi) = yi;

Một số phương pháp nội suy truyền thống được đưa ra và giải quyết rất tốt bài toán trên, điển hình là phương pháp nội suy Lagrange và

phương pháp nội suy Newton. Đa thức nội suy Lagrange rất đơn giản

4

và dễ tính, nếu các nút nội suy đã được cố định. Nhưng nếu ta bổ sung thêm nút nội suy thì quá trình tính lại phải tính lại từ đầu. Phương

pháp nội suy Newton khắc phục được nhược điểm của nội suy Lagrange ở chỗ khi thêm vào lưới nội suy một nút nội suy mới xn+1, ta chỉ cần thêm vào đa thức nội suy Pn(x) một số hạng. Tuy nhiên, khi số mốc nội suy lớn thì nội suy bằng đa thức thường xảy ra hiện tượng phù hợp

trội (overfitting) do bậc của đa thức thường tăng theo số mốc nội suy.

Hơn nữa, đa số các bài toán nội suy trong các ứng dụng thực tiễn lại là bài toán nội suy nhiều biến. Để khắc phục nhược điểm này, một phương

pháp nội suy được đề xuất bởi Powell vào năm 1987 là phương pháp nội suy hàm cơ sở bán kính (Radial basis function-RBF) có thể chuyển từ

bài toán nội suy hàm nhiều biến về nội suy hàm một biến. Hơn nữa còn cho kết quả rất tốt, đặc biệt với bài toán nội suy hàm nhiều biến trên

1.2 Một số định nghĩa và khái niệm

tập dữ liệu phân tán.

Bài toán 1.1

Cho bộ dữ liệu (xi; yi), i = 1, 2, ..., n, xi ∈ Rd; yi ∈ R, trong đó xi là các vị trí đo; yi là kết quả đo được tại vị trí xi. B1, B2, ..., Bn là các hàm cơ sở của không gian tuyến tính của các hàm liên tục d biến. Ký hiệu là:

k=1

(cid:40) n (cid:41) (cid:88) (1.1) F = span{B1, B2, ..., Bn} = ckBk; ck ∈ R

Tìm hàm Pf ∈ F sao cho

i = 1, 2, ..., n, (1.2) Pf (xi) = yi;

n (cid:88)

vì Pf ∈ F nên ta có

k=1

x ∈ Rd. (1.3) Pf (x) = ckBk(x),

Từ (1.2) và (1.3) ta có

Ac = y, (1.4)

5

trong đó  

B1(x1) B2(x1) ... Bn(x1) B2(x1) B2(x2) ... Bn(x2) A = , (1.5) ...           Bn(x1) Bn(x2) ... Bn(xn)

c = [c1, c2, ..., cn]T ; y = [y1, ..., yn]T .

Bài toán 1.1 có nghiệm duy nhất khi và chỉ khi ma trận A không suy

biến, tức là detA (cid:54)= 0. Trường hợp d = 1 (trong không gian một chiều) ta có thể chọn cơ sở

như sau:

{B1, B2, ..., Bn} = {1, x, x2, ..., xn−1}.

Tuy nhiên khi d ≥ 2 ta có kết quả sau:

Định lý 1.2.1 [10](Mairhuber - Curtis) Nếu Ω ⊂ Rd, d ≥ 2 và chứa một điểm trong thì không tồn tại không gian Haar các hàm liên tục trên Ω.

Trong đó, không gian Haar được định nghĩa như sau:

Định nghĩa 1.2.2 Cho Ω ⊂ Rd, và F ⊂ C(Ω) là không gian tuyến tính hữu hạn chiều có cơ sở là {B1, B2, ..., Bn}. Ta nói F là không gian Haar trên Ω nếu detA (cid:54)= 0 với mọi bộ tâm phân biệt {x1, x2, ..., xn} trong Ω. Trong đó ma trận

A = (Ajk)n×n; Ajk = Bk(xj); j, k = 1, 2, ..., n.

Sự tồn tại của không gian Haar đảm bảo tính khả nghịch của ma trận nội suy, nghĩa là tồn tại duy nhất nghiệm của Bài toán 1.1.

Ví dụ 1.2.3 Không gian các đa thức một biến bậc n − 1 chính là không gian Haar n chiều với tập dữ liệu (xj; yj), j = 1, ...n; xj ∈ R; yj ∈ R.

Định lí Mairhuber – Curtis cho thấy rằng nếu muốn giải được bài toán

nội suy với dữ liệu phân tán trong không gian nhiều chiều thì cơ sở cần phụ thuộc vào các vị trí dữ liệu. Để thu được các không gian xấp xỉ phụ

6

thuộc dữ liệu, chúng ta cần xét đến các hàm xác định dương và các ma trận dương.

Định nghĩa 1.2.4 (Ma trận xác định dương)

n (cid:88)

n (cid:88)

Ma trận A vuông, đối xứng, giá trị thực được gọi là xác định dương nếu dạng toàn phương tương ứng không âm:

j=1

k=1

cjckAjk ≥ 0 với c = (c1, c2, ..., cn)T ∈ Rn.

Dấu bằng chỉ xảy ra khi và chỉ khi c = (0, 0, ..., 0)T .

k=1, trong Bài toán 1.1 làm cho ma trận nội suy A

Tính chất quan trọng của ma trận xác định dương là nó có tất cả các

giá trị riêng đều dương và không suy biến. Nếu hệ cơ sở {Bk}n xác định dương thì hệ (1.4) có nghiệm duy nhất.

n (cid:88)

n (cid:88)

Định nghĩa 1.2.5 (Hàm xác định dương) Hàm liên tục Φ : Rd −→ R là xác định dương trên Rd khi và chỉ khi nó là hàm chẵn và với mọi bộ tâm phân biệt từng đôi một X = {x1, x2, ..., xn} ⊂ Rd; n ∈ N, và mọi vectơ c = (c1, c2, ...cn) ∈ Rn thì dạng toàn phương

j=1

k=1

(1.6) cjckΦ(xj − xk) ≥ 0.

Đẳng thức xảy ra khi và chỉ khi c = (0, 0, ..., 0).

Định nghĩa 1.2.6 Hàm một biến φ : [0, ∞) −→ R được gọi là xác định dương trên Rd nếu hàm nhiều biến tương ứng Φ(x) = φ(||x||) với ∀x ∈ Rd là xác định dương.

n (cid:88)

Từ định nghĩa trên và tính chất của ma trận xác định dương ta thấy có thể sử dụng các hàm xác định dương Bn = Φ(x − xk) là hệ hàm cơ sở và khi đó (1.3) trở thành

k=1

(1.7) Pf (x) = ckΦ(x − xk).

7

Ma trận nội suy A = [Ajk]n×n với Ajk = Bk(xj) = Φ(xj − xk); j, k = 1, ..., n.

Tuy nhiên việc giải bài toán nội suy trong không gian nhiều chiều là khó khăn, do đó thay vì sử dụng hàm nhiều biến Φ(x) bởi hàm một biến φ

cho tất cả số chiều d.

Định nghĩa 1.2.7 (Hàm bán kính) Hàm Φ : Rd → R được gọi là hàm bán kính nếu tồn tại hàm một biến φ : [0, ∞) → R sao cho Φ(x) = φ(||x||) với ∀x ∈ Rd.

Định nghĩa 1.2.8 (Hàm bán kính xác định dương) Cho hàm Φ : Rd → R với hàm cơ sở tương ứng là φ. Ta nói φ xác định dương trên Rd khi và chỉ khi Φ xác định dương trên Rd.

n (cid:88)

Định nghĩa 1.2.9 (Hàm bán kính xác định dương có điều kiện) Hàm chẵn, liên tục Φ : Rd → R được gọi là xác định dương có điều kiện bậc l nếu với mọi bộ tâm phân biệt từng đôi một X = {x1, x2 . . . , xn} ⊂ Rd, n ∈ N, với mọi vectơ c = {c1, c2 . . . , cn} ∈ Rn và mọi đa thức P giá trị thực bậc nhỏ hơn l, thỏa mãn

j=1

cjp(xj) = 0,

n (cid:88)

n (cid:88)

thì

j=1

k=1

cjckΦ(xj − xk) ≥ 0,

và công thức trên là đẳng thức khi và chỉ khi c là vectơ 0.

Nhận xét i) Nếu một hàm là xác định dương có điều kiện bậc l trong không gian Rd thì nó sẽ là xác định dương có điều kiện với mọi bậc lớn hơn l. Cụ thể là nếu một hàm là xác định dương (l = 0) thì sẽ là xác định dương với mọi bậc l ∈ N. ii) Ma trận A với các phần tử Aj,k = Φ(xj − xk) tương ứng với hàm chẵn, liên tục và xác định dương có điều kiện bậc l, có thể được xem

8

n (cid:88)

như là hàm xác định dương trên không gian vectơ c sao cho

j=1

cjp(xj) = 0,

trong đó p là đa thức bậc nhỏ hơn l.

Định nghĩa 1.2.10 (Hàm cơ sở bán kính (Radial basis function - RBF)) Hàm Φ : Rd → R được gọi là hàm bán kính nếu tồn tại hàm số một biến φ : [0; ∞) → R thỏa mãn:

Φ(x) = φ(r),

với r = ||x||, ||.|| là một chuẩn nào đó trong Rd (ta thường dùng chuẩn Ơcơlit). φ được gọi là hàm cơ sở bán kính.

Trong khuôn khổ luận văn này tôi trình bày một số hàm cơ sở bán

kính thông dụng, với r = (cid:107)x − xk(cid:107) .

Tên hàm Tên viết tắt Định nghĩa √ Multiquadric MQ φmq(r) =

√ Inverse multiquadric IMQ φimq(r) = 1 + r2 1 1 + r2

Gaussian Gauss

Cauchy Cauchy φc(r) = φg(r) = e−r2 1 1 + r2

Bảng 1.1: Bảng một số hàm cơ sở bán kính dùng trong luận văn.

Hàm cơ sở bán kính φ(x) là xác định dương nếu ta nhân r với một số

dương ε thì φ(rε) vẫn là xác định dương. Khi các mốc nội suy xác định thì giải pháp tối ưu là đưa vào hàm Φk một tham số hình dạng εk. Như vậy ta cần tìm εk để bài toán thỏa mãn điều kiện nội suy, đồng thời chất lượng nội suy là tốt nhất. Khi đó εk còn gọi là tham số tỉ lệ (scaling) của hàm cơ sở bán kính vì nó dùng để điều chỉnh độ rộng của miền ảnh hưởng của hàm cơ sở φ. Khi ||x − xk|| > ε − εk thì giá trị hàm Φk(x) là rất nhỏ không có ý nghĩa vì nó gần triệt tiêu. Vì vậy ta nói hàm bán kính này chỉ có ảnh hưởng địa

9

phương. * Tham số hình dạng cho một số hàm cơ sở bán kính [10].

Tên hàm Tên viết tắt Biểu thức tham số hóa hình dạng √ Multiquadric MQ φmq(εr) =

√ Inverse mulhquadric IMQ φimq(εr) = ε2 + r2 1 ε2 + r2

Gaussian Gauss φg(εr) = e−(εr)2

Cauchy Cauchy φc(εr) = 1 ε2 + r2

Bảng 1.2: Bảng một số hàm cơ sở bán kính với tham số hình dạng ε > 0.

Cho miền Ω trong không gian Euclide Rd với biên ∂Ω. Ta sẽ dùng

thuật ngữ "Tâm" như là một điểm thuộc miền Ω.

n (cid:88)

Định nghĩa 1.2.11 (Véc tơ trọng số (Stencil)) Cho D là toán tử vi phân tuyến tính và X = {x1, x2 . . . , xn} là bộ tâm phân tán đã được chọn trong không gian Rd. Một xấp xỉ của toán tử D,

i=1

Du(x) ≈ (1.8) wi(x)u(xi),

được xác định bởi các trọng số wi = wi(x). Khi đó w = [w1, w2 . . . , wn]T được gọi là véctơ trọng số hay còn được gọi là stencil đối với toán tử vi phân D.

Định nghĩa 1.2.12 (Bộ tâm nằm trong miền Ξint) Ξint là bộ tâm nằm trong miền được định nghĩa bởi Ξint = Ξ \ ∂Ξ, trong đó, Ξ là các tâm trên toàn bộ miền và ∂Ξ là các tâm nằm trên biên.

Định nghĩa 1.2.13 (Bộ tâm hỗ trợ cho véctơ trọng số Ξζ) Với mỗi ζ ∈ Ξint, ta chọn bộ tâm bao gồm ζ và các tâm nằm trong lân cận địa phương của nó. Bộ tâm này được gọi là bộ tâm để tính véctơ trọng số và được kí hiệu là Ξζ.

10

1.3 Lược đồ phương pháp sai phân hữu hạn giải phương trình đạo hàm riêng trong không gian 3 chiều

Bài toán truyền nhiệt dừng trong không gian 3 chiều

Cho các số a, b, c, d, e, f với a < b, c < d, e < f. Xét trong không gian với hệ tọa độ vuông góc Oxyz một khối hộp chữ nhật Ω có các mặt song

song với các mặt phẳng tọa độ (Oxy), (Oyz), (Ozx).

Ω = {(x, y, z)|a < x < b, c < y < d, e < z < f }

có mặt biên khép kín và ký hiệu là ∂Ω. Xét bài toán: Tìm hàm số u(x, y, z) thỏa mãn phương trình Poission:

(1.9) ∆u ≡ ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 = f (x, y, z), (x, y, z) ∈ Ω

và điều kiện biên

u(x, y, z) = g(x, y, z), (x, y, z) ∈ ∂Ω (1.10)

trong đó f (x, y, z) và g(x, y, z) là các hàm số cho trước.

Giả sử bài toán (1.9) – (1.10) có nghiệm duy nhất u = u(x, y, z) đủ

N là bước đi là bước đi theo z. Đặt

trơn trong Ω = Ω ∪ ∂Ω.

P

Bước 1: Xây dựng lưới sai phân và hàm lưới *) Chọn các số nguyên N > 1, M > 1, P > 1, đặt h = b−a M gọi là bước đi theo y, p = f −e

theo x, k = d−c xi = a + ih, yj = c + jk, zn = e + np.

Mỗi điểm (xi, yj, zn) gọi là một nút lưới, còn được ký hiệu là (i, j, n). Tập Ωhkp = {(xi, yj, zn)|1 ≤ i ≤ N − 1, 1 ≤ j ≤ M − 1, 1 ≤ n ≤ P − 1} gọi là tập các nút trong. Tập ∂Ωhkp = {(x0, yj, zn), (xN , yj, zn), (xi, y0, zn), (xi, yM , zn), (xi, yj, z0), . . . (xi, yj, zP )} gọi là tập các nút biên. Tập Ωhkp = Ωhkp ∪ ∂Ωhkp gọi là một lưới sai phân trên Ω.

*) Hàm số v xác định tại mọi nút của lưới Ωhkp gọi là một hàm lưới

trên Ωhkp. Giá trị của hàm lưới v tại nút (i, j, n) viết là vn ij.

*) Đạo hàm cấp 1 tiến, đạo hàm cấp 1 lùi của hàm lưới v xác định

trên Ωhkp:

i−1j

ij =

ij =

ij−1

vx xác định bởi (vx)n ; vx xác định bởi (vx)n

ij =

ij =

ij

vy xác định bởi (vy)n ; vy xác định bởi (vy)n

ij =

ij =

11 i+1j − vn vn ij h ij+1 − vn vn ij k vn+1 ij − vn ij p

ij − vn vn h ij − vn vn k ij − vn−1 vn p

vz xác định bởi (vz)n ; vz xác định bởi (vz)n

*) Đạo hàm cấp 2 của hàm lưới v xác định trên Ωhkp:

i+1j − (vx)n ij

ij = ((vx)x)n

ij =

i+1j − 2vn

ij + vn

i−1j]

(vx)n = (vxx)n h

ij+1 − (vy)n ij

ij = ((vy)y)n

ij =

ij+1 − 2vn

ij + vn

ij−1]

(vy)n = (vyy)n

ij − (vz)n ij

ij = ((vz)z)n

ij =

ij + vn−1

ij − 2vn

ij

k (vz)n+1 = ] (vzz)n p 1 h2 [vn 1 k2 [vn 1 p2 [vn+1

Ta sẽ tính gần đúng giá trị của nghiệm u(x, y, z) tại các nút (xi, yj, zn) và kí hiệu giá trị gần đúng đó là vn ij:

vn ij ≈ u(xi, yj, zn)

Bước 2: Xây dựng bài toán sai phân

Áp dụng công thức khai triển Taylor cho:

u(xi+1, yj, zn) = u(xi + h, yj, zn)

+ = u(xi, yj, zn) + h h2 2! ∂2u ∂x2 + h3 3! ∂3u ∂x3 + 0(h4)

∂u ∂x u(xi−1, yj, zn) = u(xi − h, yj, zn)

+ = u(xi, yj, zn) − h ∂u ∂x h2 2! ∂2u ∂x2 − h3 3! ∂3u ∂x3 + 0(h4)

Vậy có:

= u(xi+1, yj, zn) − 2u(xi, yj, zn) + u(xi−1, yj, zn) h2 ∂2u ∂x2 + 0(h2)

Một cách tương tự:

= u(xi, yj+1, zn) − 2u(xi, yj, zn) + u(xi, yj−1, zn) k2

= u(xi, yj, zn+1) − 2u(xi, yj, zn) + u(xi, yj, zn−1) p2 ∂2u ∂y2 + 0(k2) ∂2u ∂z2 + 0(p2)

12

Vậy có:

+

+ + u(xi+1, yj, zn) − 2u(xi, yj, zn) + u(xi−1, yj, zn) h2 u(xi, yj+1, zn) − 2u(xi, yj, zn) + u(xi, yj−1, zn) k2

= ∆u + 0(h2 + k2 + p2) u(xi, yj, zn+1) − 2u(xi, yj, zn) + u(xi, yj, zn−1) p2

(1.11)

i−1j

ij + vn

ij−1

ij

Gọi v là một hàm lưới xác định tại mọi nút của lưới Ωhkp. Ta đặt:

i+1j − 2vn vn h2

ij + vn ij+1 − 2vn vn k2

ij + vn−1 vn+1 ij − 2vn p2

+ + ∆hkp ≡

(1.12)

Khi đó (1.11) chứng tỏ:

(1.13) ∆hkpu = ∆u + 0(h2 + k2 + p2)

Số hạng 0(h2 + k2 + p2) là một vô cùng bé. Ta nói toán tử ∆hkp xấp xỉ toán tử ∆. Do đó khiến ta thay phương trình vi phân (1.9) - (1.10) bằng

phương trình sai phân sau:

ij, f n

ij = f (xi, yj, zn), (xi, yj, zn) ∈ Ωhkp

ij + vn

i−1j

ij−1

ij

(1.14) ∆hkpv = f n

ij + vn ij+1 − 2vn vn k2

+ + Tức là: i+1j − 2vn vn h2 vn+1 ij + vn−1 ij − 2vn p2

(1.15) = f (xi, yj, zn), (xi, yj, zn) ∈ Ωhkp

và thay điều kiện biên (1.10) bằng điều kiện:

(1.16) vn ij = g(xi, yj, zn), (xi, yj, zn) ∈ ∂Ωhkp

chúng ta được bài toán sai phân hoàn chỉnh: Tìm hàm lưới v = v(xi, yj, zn) thỏa mãn phương trình sai phân (1.15) – (1.16).

Bước 3: Giải bài toán sai phân (1.15) – (1.16) Có rất nhiều phương pháp để giải hệ phương trình đại số tuyến tính

(1.15) – (1.16), ví dụ như phương pháp Gauss, phương pháp luân phương

13

ij−1, vn+1

i+1j, vn

ij+1, vn

i−1j, vn

ij

,

ẩn, phương pháp nhẩy ô, phương pháp một chiều địa phương, phương pháp sai phân hiện, phương pháp sai phân ẩn,.... , vn−1 Trong hệ phương trình (1.15) có tới 7 ẩn vn ij vn ij nên việc giải nghiệm đúng là vô cùng phức tạp, do đó ta chỉ có thể chọn phương pháp giải nghiệm gần đúng, tức là tìm được các nghiệm xấp xỉ vn ij tại các nút trong (i, j, n) của nghiệm chính xác u = u(x, y, z).

14

1.4 Kết luận

Trong chương này, chúng tôi đã trình bày lại một cách chi tiết và có

hệ thống các kiến thức liên qua đến luận văn: Cơ sở của bài toán nội suy hàm số với dữ liệu phân tán; nội suy với dữ liệu phân tán; khái

niệm về ma trận xác định dương; hàm xác định dương; hàm bán kính; hàm bán kính xác định dương; hàm bán kính xác định dương có điều

kiện; khái niệm về hàm cơ sở bán kính (RBF); khái niệm về véctơ trọng số(stencil); khái niệm về bộ tâm nằm trong miền Ξint; khái niệm về bộ tâm hỗ trợ cho véctơ trọng số Ξζ. Trình bày lại lược đồ phương pháp sai phân hữu hạn giải phương trình đạo hàm riêng trong không gian ba

chiều.

15

Chương 2

Phương pháp không lưới RBF –

FD giải phương trình Poisson

trong không gian ba chiều

Trong chương này tôi trình bày một số phương pháp xác định véctơ trọng số: Véctơ trọng số từ vi phân số trên các tâm phân bố không đều;

véctơ trọng số RBF dựa trên nội suy hàm cơ sở bán kính, trình bày lại một cách chi tiết kết quả của bài báo [15] về Octant - sự lựa chọn

véctơ trọng số cho phương pháp sai phân không lưới giải phương trình Poisson trong không gian ba chiều; Và đưa ra lược đồ RBF – FD giải

bài toán Elliptic trong không gian ba chiều. Nội dung chính của chương

2.1 Bài toán mở đầu

này được tham khảo chủ yếu từ các tài liệu [1, 6, 7, 10], và từ các bài báo [14, 15].

Xét bài toán Dirichlet với phương trình Poisson:

Tìm hàm u : Ω → R sao cho

∆u = f trong Ω, (2.1) u = g trên ∂Ω.

trong đó ∆ là toán tử Laplace, Ω ⊂ R3 là miền hình học, f là hàm giá trị thực được xác định trong miền Ω, và g là một hàm giá trị thực được

16

xác định trên ∂Ω.

Phương pháp không lưới RBF – FD được trình bày như sau

Cho Ξ là một tập các nút rời rạc, trong đó Ξ ⊂ Ω và Ξint := Ξ ∩ Ω, ∂Ξ := Ξ ∩ ∂Ω. Đối với mỗi ζ ∈ Ξint, chọn Ξζ := {ξ0, ξ1, ξ2, . . . , ξk} ⊂ Ξ với ξ0 = ζ, và xấp xỉ ∆u (ζ) bởi công thức vi phân số

ξ∈Ξζ

(cid:88) ∆u (ζ) ≈ (2.2) wζ,ξu (ξ) , ζ ∈ Ξint.

Khi đó, bài toán (2.1) có thể được viết lại như sau

wζ,ξ ˆu (ξ) = f (ζ) , ζ ∈ Ξint; (cid:80) ξ∈Ξζ (2.3)

ˆu (ξ) = g (ξ) , ξ ∈ ∂Ξ,

2.2 Véctơ trọng số từ nội suy hàm cơ sở bán kính

trong đó ˆu(ξ) là xấp xỉ của nghiệm u đối với bài toán (2.1) tại các điểm ξ ∈ Ξ. Giải hệ phương trình tuyến tính (2.3), chúng ta thu được vectơ [ˆu(ξ)]ξ∈Ξ. Một phương pháp RBF – FD bất kỳ được định nghĩa bởi thuật toán chọn tâm, thuật toán này chọn các tập hỗ trợ Ξζ ⊂ Ξ, ζ ∈ Ξint và cách tính trọng số vi phân số (stencil) wζ,ξ, ζ ∈ Ξint, ξ ∈ Ξζ. Các trọng số này có thể tìm được nhờ công thức (2.2).

bố không đều

2.2.1 Véctơ trọng số từ vi phân số trên các tâm phân

m (cid:88)

Cho bộ dữ liệu phân tán X = {x1, x2 . . . , xn} ⊂ Rd với giá trị hàm tương ứng u|X = (cid:2)u(x1), u(x2), . . . , u(xn)(cid:3)T . Giả sử s(x) là một xấp xỉ của hàm u(x) từ tập dữ liệu X, dưới dạng

i=1 trong đó si(x), i = 1, 2, . . . , m (m ∈ Z, m ≥ n và m hữu hạn) là các hàm cơ sở của không gian tuyến tính P nào đó và véctơ hệ số

s(x) = (2.4) aisi(x),

17

j=1

a = [a1, a2, . . . , am]T phụ thuộc tuyến tính vào u|X, n (cid:88) i = 1, 2, . . . , m. (2.5) ai = biju(xj),

i=1,j=1, khi đó ta có a = B · u|X .

Ký hiệu B = [bij]m,n

m (cid:88)

Vì D là toán tử vi phân tuyến tính nên ta có,

Du(x) ≈ Ds(x) = aiDsi(x)

i=1 n (cid:88)

i=1

j=1 n (cid:88)

(cid:17) (cid:16) m (cid:88) = bijDsi(x) u(xj)

j=1

= wju(xj),

m (cid:88)

và ta nhận được công thức vi phân số (1.8) với véctơ trọng số w = [w1, w2, . . . , wn]T được xác định bởi

i=1

j = 1, 2, . . . , n, wj = bijDsi(x),

hoặc dưới dạng ma trận

i=1.

(2.6) w = BT · [Dsi(x)]m

Trường hợp riêng, khi m = n và các hệ số aj thu được bằng cách giải

n (cid:88)

bài toán nội suy không suy biến

j=1

i,j=1. Vì vậy, véctơ

i = 1, 2, . . . , n. s(xi) = ajsj(xi) = u(xi),

Khi đó rõ ràng B = [S|X]−1, trong đó S|X := [sj(xi)]n trọng số w được xác định bởi

i=1.

w = [S|X]−T · [Dsi(x)]n

Và có thể tìm bằng cách giải hệ phương trình

i=1,

(2.7) [S|X]T w = [Dsi(x)]n

18

n (cid:88)

tức là,

i=1

j = 1, 2, . . . , n. wisj(xi) = Dsj(x),

Điều quan trọng cần lưu ý là từ các công thức trên ta thấy rằng véctơ

trọng số được xác định duy nhất bởi toán tử vi phân D và các hàm cơ sở s1, s2 . . . , sn, do đó nó không thay đổi nếu thực hiện phép biến đổi tuyến tính không suy biến các hàm cơ sở.

Rõ ràng có nhiều phương pháp xấp xỉ khác nhau, ví dụ như phương

pháp bình phương tối thiểu hoặc tựa nội suy kiểu (2.4)–(2.5) có thể dẫn đến các công thức tính véctơ trọng số dưới dạng (2.6).

m (cid:88)

n (cid:88)

i=1

j=1

Vì (cid:16) (cid:17) Du(x) − u(x) − , wju(xj) = D aisi(x)

nên độ chính xác của các véctơ trọng số liên quan trực tiếp đến độ chính

xác của phương pháp xấp xỉ.

2.2.2 Nội suy không có thành phần đa thức

n (cid:88)

Cho bộ tâm phân biệt từng đôi một X = {x1, x2, . . . , xn} ⊂ Rd, u : Rd → R là hàm liên tục và đủ trơn. Giả sử φ : R+ → R là hàm xác định dương và đủ trơn. Khi đó hàm nội suy cơ sở bán kính [6, 7, 10] s(x) của hàm u(x) được viết dưới dạng

j=1

s(x) = Φ(x) := φ((cid:107)x(cid:107)), (2.8) ajΦ(x − xj),

i = 1, . . . , n, (2.9) s(xi) = u(xi),

trong đó aj được được chọn sao cho thỏa mãn điều kiện nội suy (2.9).

n (cid:88)

Từ (2.8) và (2.9), ta có:

j=1

i = 1, . . . , n. (2.10) ajΦ(xi − xj) = u(xi),

Kí hiệu

i,j=1, u|X = [u(x1), u(x2), . . . , u(xn)]T , a = [a1, a2, . . . , an]T .

Φ|X = [Φ(xi − xj)]n,n

19

Khi đó, ta có thể viết (2.10) dưới dạng ma trận

Φ|X a = u|X .

Vì φ là hàm xác định dương nên ma trận Φ|X là xác định dương với bộ tâm X phân biệt từng đôi một. Do đó, véctơ a được xác định duy nhất bởi

(2.11) a = (cid:2) Φ|X (cid:3)−1 u|X .

Hàm nội suy cơ sở hàm bán kính s(x) là một xấp xỉ tốt của hàm u(x) nếu hàm u(x) đủ trơn và các tâm x1, x2, . . . , xn ∈ Rd đủ dầy trong lân cận của x. Hơn nữa, đạo hàm của hàm s(x) cũng là xấp xỉ tốt với đạo

hàm của hàm u(x) nếu hàm φ đủ trơn. Vì vậy, nếu D là toán tử vi phân

n (cid:88)

tuyến tính thì Ds(x) là xấp xỉ của Du(x) được xét dưới dạng

j=1

Du(x) ≈ Ds(x) = ajDΦ(x − xj)

(2.12) = aT D Φ(x − ·)|X (cid:2) Φ|X = u|T X (cid:3)−1D Φ(x − ·)|X .

(cid:3)T . Ta đặt Ký hiệu w = (cid:2)w1, w2, . . . , wn

(2.13) w = (cid:2) Φ|X (cid:3)−1D Φ(x − ·)|X ,

trong đó D Φ(x − ·)|X = (DΦ(x − x1), . . . , DΦ(x − xn))T .

n (cid:88)

Từ (2.12) và (2.13) ta có thể viết

i=1

Du(x) ≈ Ds(x) = (2.14) wiu(xi),

n (cid:88)

ta nhận được công thức vi phân số (1.8) với véctơ trọng số w = [w1, w2, . . . , wn]T được xác định bởi (2.13). Quan sát công thức (2.13) ta thấy w là nghiệm của hệ phương trình

j=1

i = 1, 2, . . . , n, (2.15) wjΦ(xi − xj) = DΦ(x − xi),

Điều này có nghĩa là véctơ trọng số w được cho bởi các hệ số của nội suy hàm cơ sở bán kính với dữ liệu cho bởi hàm D Φ(x − ·)|X. Vì

20

Φ(xi − xj) = Φ((x − xi) − (x − xj)) nên nếu ta nội suy hàm DΦ tại các tâm x − xj, j = 1, 2, . . . , n thì ta thu được các hệ số nội suy như trong công thức (2.13) và đó chính là các véctơ trọng số. Vì vậy chúng ta có phương pháp tính véctơ trọng số như sau:

Mệnh đề 2.2.1 Cho bộ tâm phân biệt từng đôi một X = {x1, x2, . . . , xn} ⊂ Rd, u : Rd → R là hàm liên tục và đủ trơn. Giả sử φ : R+ → R là hàm xác định dương và đủ trơn, D là toán tử vi phân tuyến tính và hàm nội suy cơ sở bán kính s(x) của hàm u(x) được biểu diễn dưới dạng (2.8)–(2.9). Khi đó, véctơ trọng số w của vi phân số tại

x được tìm bằng cách giải hệ phương trình (2.15), hay véctơ trọng số là các hệ số của nội suy hàm cơ sở bán kính với dữ liệu được cho bởi hàm D Φ(x − ·)|X.

2.2.3 Nội suy có thành phần đa thức

Cho bộ tâm phân biệt từng đôi một X = {x1, x2, . . . , xn} ⊂ Rd, u : Rd → R là hàm liên tục và đủ trơn. Giả sử φ : R+ → R là hàm xác định dương hoặc hàm xác định dương có điều kiện bậc cao nhất là (cid:96) + 1

và đủ trơn.

Khi đó s(x) là hàm nội suy cơ sở bán kính của hàm u(x) được viết

n (cid:88)

dưới dạng

j=1

s(x) = Φ(x) := φ((cid:107)x(cid:107)), (2.16) ajΦ(x − xj) + ℘(x),

i = 1, 2, . . . , n, ℘ là đa thức bậc (cid:96). (2.17) s(xi) = u(xi),

n (cid:88)

Từ (2.16) và (2.17), ta có thể viết

j=1

i = 1, 2, . . . , n. (2.18) ajΦ(xi − xj) + ℘(xi) = u(xi),

Vì φ là hàm xác định dương hoặc xác định dương có điều kiện bậc cao nhất là (cid:96) + 1 nên các hệ số aj và đa thức ℘(x) được xác định duy nhất

21

n (cid:88)

[10, Định lý 8.21] bởi các điều kiện

j=1 n (cid:88)

i = 1, 2, . . . , n, (2.19) ajΦ(xi − xj) + ℘(xi) = u(xi),

j=1

= 0, p là đa thức bậc (cid:96). (2.20) ajp(xj)

k (cid:88)

Giả sử rằng p1, p2, . . . , pk là cơ sở của các đa thức bậc (cid:96), ta ký hiệu

i,j=1,

℘(x) = a = [a1, a2, . . . , an]T , cjpj(x), Φ|X = [Φ(xi − xj)]n,n

j=1 c = [c1, c2, . . . , ck]T .

i=1,j=1,

P |X = [pj(xi)]n,k u|X = [u(x1), . . . , u(xn)]T ,

Khi đó, chúng ta có thể viết hệ (2.19)–(2.20) dưới dạng ma trận

X a

Φ|X a + P |X c = u|X , P |T = 0.

Do đó, véctơ a và c được xác định duy nhất bởi

(cid:34) (cid:35) (cid:34) (cid:35)−1 (cid:34) (cid:35)

= a c u|X 0 Φ|X P |X P |T 0 X

Quay lại Phần 2.2.1, chúng ta thấy rằng ở đây m = n + k, si = φi, i = 1, 2, . . . , n, sn+j = pj, j = 1, 2, . . . , k, và ma trận B thu được từ ma (cid:35)−1 (cid:34)

bằng cách bỏ đi k cột cuối cùng. trận ΦX P |X P |T 0 X

Hàm nội suy cơ sở bán kính s(x) là một xấp xỉ tốt của hàm u(x) nếu hàm u(x) đủ trơn và các điểm x1, x2, . . . , xn ∈ Rd đủ dầy trong lân cận của x. Hơn nữa, đạo hàm của hàm s(x) cũng là xấp xỉ tốt với đạo hàm của hàm u(x) nếu hàm φ đủ trơn. Những khẳng định này có thể xem

trong. Vì vậy, một xấp xỉ của Du(x), trong đó D là toán tử vi phân

22

k (cid:88)

n (cid:88)

tuyến tính có thể được xét dưới dạng,

j=1 (cid:35)

j=1 (cid:34)

Du(x) ≈ Ds(x) = ajDΦ(x − xj) + cjDpj(x)

(cid:35)T (cid:34)

= DΦ(x − ·)|X DP (x) a c

(cid:34) (cid:35)T (cid:34) (cid:35)−1 (cid:34) (cid:35)

= u|X 0 DΦ(x − ·)|X DP (x)

n (cid:88)

(cid:34) (cid:35)T (cid:34) Φ|X P |X P |T 0 X (cid:35)

i=1

= = (2.21) wiu(xi), w v u|X 0

trong đó, DP (x) = (Dp1(x), Dp2(x), . . . , Dpk(x))T và

(cid:35) (cid:35)−1 (cid:34) (cid:34) (cid:35) (cid:34) w (2.22) = DΦ(x − ·)|X DP (x) v Φ|X P |X P |T 0 X

Bằng cách nhân cả hai vế của (2.22) với

(cid:34) (cid:35)

, Φ|X P |X P |T 0 X

dẫn đến véctơ trọng số w thỏa mãn

X w

Φ|X w + P |X v = DΦ(x − ·)|X, P |T = DP (x)

n (cid:88)

hay véctơ trọng số w được tìm từ hệ phương trình

j=1 n (cid:88)

(2.23) wjΦ(xi − xj) + ℘(xi) = DΦ(x − xi), i = 1, n,

j=1

= Dp(x) (2.24) wjp(xj)

Những lập luận trên có thể dẫn đến

Mệnh đề 2.2.2 Cho bộ tâm phân biệt từng đôi một X = {x1, x2, . . . , xn} ⊂ Rd, u : Rd → R là hàm liên tục và đủ trơn. Giả sử φ : R+ → R là hàm xác định dương hoặc xác định dương có điều kiện

23

bậc cao nhất là (cid:96) + 1 và đủ trơn, D là toán tử vi phân tuyến tính và hàm nội suy cơ sở bán kính s(x) của hàm u(x) được biểu diễn dưới dạng

(2.16)–(2.17). Khi đó, véctơ trọng số w của vi phân số tại x được tìm từ hệ phương trình (2.23)–(2.24).

Nếu (cid:96) ≤ 1 có nghĩa là φ là hàm xác định dương hoặc xác định dương có

điều kiện bậc cao nhất bằng 1 và p là đa thức bậc 0. Trong trường hợp này, bài toán (2.23)–(2.24) giống bài toán nội suy (2.19)–(2.20) và bài

n (cid:88)

toán (2.23)–(2.24) được viết lại như sau

j=1 n (cid:88)

(2.25) wjΦ(xi − xj) + v = DΦ(x − xi), i = 1, n,

j=1

= 0 (2.26) wj

Nhận xét này dẫn đến phát biểu sau

Mệnh đề 2.2.3 Cho bộ tâm phân biệt từng đôi một X = {x1, x2, . . . , xn} ⊂ Rd, u : Rd → R là hàm liên tục và đủ trơn. Giả sử φ là hàm xác định dương hoặc xác định dương có điều kiện bậc cao nhất là 1 và đủ trơn, D là toán tử vi phân tuyến tính.

Khi đó véctơ trọng số w của vi phân số tại x được cho bởi các hệ số của nội suy hàm cơ sở bán kính φ với các tâm x − xj, j = 1, 2, . . . , n và dữ liệu được cho bởi hàm DΦ(x − xi), i = 1, 2, . . . , n. Cụ thể là w là nghiệm của hệ phương trình tuyến tính (2.25)–(2.26).

2.2.4 Véctơ trọng số RBF trong không gian ba chiều

Cho bộ tâm phân biệt từng đôi một X = {x1, x2, ..., xn} ⊂ R3,

n (cid:88)

u : R3 → R là hàm liên tục và đủ trơn. Giả sử φ : R+ → R là hàm xác định dương và đủ trơn [6]. Khi đó hàm nội suy RBF s(x) của hàm u(x) được xác định bởi công thức:

j=1

s(x) = (2.27) ajΦ(x − xj), Φ(x) := φ(||x||),

24

sao cho

(2.28) s(xi) = u(xi), i = 1, 2, ..., n.

n (cid:88)

Các hằng số aj được chọn để điều kiện nội suy (2.28) được thoả mãn. Từ (2.27) − (2.28) ta có:

j=1

i = 1, 2, ..., n (2.29) ajΦ(xi − xj) = u(xi),

i,j=1,

Ký hiệu:

Φ|X = [Φ(xi − xj)]n,n u|X = [u(x1), u(x2), ..., u(xn)]T , a = [a1, a2, ..., an]T ,

khi đó ta có thể viết (2.29) dưới dạng ma trận:

Φ|Xa = u|X.

Vì φ là hàm xác định dương nên ma trận Φ|X là xác định dương với bộ tâm X phân biệt từng đôi một. Do đó, véctơ a được xác định duy

nhất bởi:

(2.30) a = [Φ|X]−1u|X.

Thay vì tính đạo hàm của hàm u(x), ta tính đạo hàm của hàm s(x). Để

làm được việc này, ta sẽ áp toán tử đạo hàm D lên hàm s(x) được biểu diễn bởi công thức (2.27) và ta có:

n (cid:80) j=1

Ds(x) = ajDΦ(x − xj)

(2.31) = aT DΦ(x − .)|X

X[Φ|X]−1DΦ(x − .)|X.

= u|T

Ký hiệu: w = [w1, w2, ..., wn]T . Ta đặt

(2.32) w = [Φ|X]−1DΦ(x − .)|X,

25

trong đó,

DΦ(x − .)|X = (DΦ(x − x1), ..., DΦ(x − xn))T .

n (cid:88)

Từ (2.31) − (2.32) ta có thể viết:

i=1

Ds(x) = (2.33) wi(x)u(xi).

n (cid:88)

Ta nhận được công thức đạo hàm số (2.33) với véctơ trọng số w = [w1, w2, ..., wn]T được xác định bởi (2.32). Quan sát công thức (2.32) ta có thể thấy véctơ trọng số w là nghiệm của phương trình:

j=1

(2.34) wjΦ(xi − xj) = DΦ(x − xi), i = 1, 2, ..., n.

Điều này có nghĩa là véctơ trọng số w được cho bởi các hệ số của nội suy hàm cơ sở bán kính với dữ liệu cho bởi hàm DΦ(x − .)|X. Vì Φ(xi − xj) = Φ((x − xi) − (x − xj)) nên nếu ta nội suy hàm DΦ tại các tâm x − xj, j = 1, 2, ..., n, thì ta thu được các hệ số nội suy như trong công thức (2.32) và đó chính là véctơ trọng số. Vì vậy chúng ta có

phương pháp tính véctơ trọng số như sau:

Mệnh đề 2.2.4 Cho bộ tâm phân biệt từng đôi một X = {x1, x2, ..., xn} ⊂ R3, u : R+ → R là hàm liên tục và đủ trơn, D là toán tử đạo hàm và hàm nội suy cơ sở bán kính s(x) của hàm u(x) được biểu diễn dưới dạng (2.27) − (2.28). Khi đó véctơ trọng số w của

đạo hàm tại x được tìm bằng cách giải hệ phương trình (2.34), hay véctơ trọng số là các hệ số của nội suy hàm cơ sở bán kính với dữ liệu được cho bởi hàm DΦ(x − .)|X.

Giả sử cho bộ tâm Ξ ⊂ Ω, để tính được đạo hàm tại điểm ζ ∈ Ξ tại theo

công thức (2.33), ta chọn bộ tâm cho nội suy hàm RBF, từ bây giờ ta sẽ ký hiệu bộ tâm này là Ξζ và gọi là bộ tâm hỗ trợ tính véctơ trọng số hay hỗ trợ tính đạo hàm, Ξζ bao gồm ζ và một số điểm nằm trong vị trí lân cận với nó.

26

2.3 Thuật toán chọn tâm dựa trên các góc khối

ζ := (cid:11).

Thuật toán 2.3.1 (Thuật toán 8 góc khối) Đầu vào: Ξ, ζ ∈ Ξint. Đầu ra: Ξζ, Ξ(cid:48) ζ. Tham số: m ≥ 16 (số điểm được chọn ban đầu gồm cả ζ). Khởi tạo: Ξζ := {ζ} , Ξ(cid:48)

I. Tìm m điểm ξ1, ..., ξm thuộc Ξ \ {ζ} và gần ζ nhất.

II. Phân hoạch các điểm ξ1, ..., ξm vào 8 góc khối Oj = {ξj1, ξj2, ...} ,

j = 1, 2, ..., 8 tương ứng với 8 góc khối, thỏa mãn (cid:107) ξj1 − ζ (cid:107)≤(cid:107) ξj2 − ζ (cid:107)≤ · · · ≤(cid:107) ξj8 − ζ (cid:107) .

III. Cho j = 1 đến 8

a. Nếu #Oj = 1 thì Ξζ := Ξζ ∪ {ξj1} .

b. Ngược lại, nếu #Oj > 1 thì Ξζ := Ξζ ∪ {ξj1, ξj2} .

ζ := Ξ(cid:48)

ζ

(cid:83) (cid:8)ξ(cid:48)(cid:9),

IV. Với mỗi ξ ∈ Ξ\{ζ}, xét đoạn thẳng (ζ, ξ) = {ζ + α(ξ − ζ) : 0 < α < 1} . V. Nếu (ζ, ξ) (cid:84) ∂Ω (cid:54)= ∅ thì Ξζ := Ξζ \ {ξ} (cid:83) (cid:8)ξ(cid:48)(cid:9) và Ξ(cid:48) trong đó ξ(cid:48) là điểm thuộc (ζ, ξ) (cid:84) ∂Ω và gần ζ nhất.

Sau khi áp dụng Thuật toán 2.3.1 cho tất cả các điểm ζ ∈ Ξint, ta cập nhật Ξ bởi công thức

(cid:48) Ξ ζ

ζ∈Ξint

(cid:88) (2.35) Ξ := Ξ ∪

Thuật toán 2.3.2 (Thuật toán 16 góc khối) Đầu vào: Ξ, ζ ∈ Ξint. Đầu ra: Ξζ, Ξ(cid:48) ζ. Tham số: m ≥ 16 (số điểm được chọn lựa ban đầu gồm cả ζ). Khởi tạo: Ξζ := {ζ} , Ξ(cid:48) ζ := (cid:11). I. Tìm m điểm ξ1, ..., ξm thuộc Ξ \ {ζ} và gần ζ nhất. II. Phân hoạch các điểm ξ1, ..., ξm vào 16 tập Oj = {ξj1, ξj2, ...} , j = 1, 2, ..., 16 tương ứng với 16 góc khối, thỏa mãn (cid:107) ξj1 − ζ (cid:107)≤(cid:107) ξj2 − ζ (cid:107)≤ · · · ≤(cid:107) ξj16 − ζ (cid:107) . III. Cho j = 1 đến 16:

27

ζ

ζ := Ξ(cid:48)

Nếu #Oj ≥ 1 thì Ξζ := Ξζ ∪ {ξj1} . IV. Với mỗi ξ ∈ Ξ \ {ζ}, xét đoạn thẳng (ζ, ξ). (cid:83) (cid:8)ξ(cid:48)(cid:9), V. Nếu (ζ, ξ) (cid:84) ∂Ω (cid:54)= ∅ thì Ξζ := Ξζ \ {ξ} (cid:83) (cid:8)ξ(cid:48)(cid:9) và Ξ(cid:48) trong đó ξ(cid:48) là điểm thuộc (ζ, ξ) (cid:84) ∂Ω và gần ζ nhất. Tương tự phía trên, sau khi áp dụng Thuật toán 2.3.2 cho tất cả các điểm ζ ∈ Ξint, ta cập nhật Ξ bởi công thức (2.35).

Tương tự phía trên, sau khi áp dụng Thuật toán 2.3.2 cho tất cả

2.4

Lược đồ RBF-FD giải bài toán Elliptic trong

không gian ba chiều

các điểm ζ ∈ Ξint, ta cập nhật Ξ bởi công thức (2.35).

Đầu vào: Ξ ⊂ Ω là tập các tâm rời rạc, Ξint := Ξ ∩ Ω là tập các tâm

trong miền và ∂Ξ := Ξ ∩ ∂Ω là tập các tâm nằm trên biên.

Đầu ra: Nghiệm xấp xỉ ˆu(ξ) của hệ phương trình (2.1) .

Bước I. Đối với mỗi ζ ∈ Ξint:

1. Chọn tập hộ trợ Ξζ bởi Thuật toán 2.3.1 hoặc Thuật

toán 2.3.2 phía trên.

2. Tính véctơ trọng số wζ,ξ với ξ ∈ Ξζ bằng cách giải hệ

phương trình đại số tuyến tính (2.34).

3. Thay véctơ wζ,ξ vào hệ phương trình (2.3).

Bước II. Giải hệ phương trình (2.3), ta thu được nghiệm xấp xỉ ˆu(ξ).

2.5 Thử nghiệm số

Bước III. Đánh giá sai số của nghiệm xấp xỉ ˆu(ξ) vừa tìm được.

Chúng tôi xem xét hai bài toán kiểu (2.1). Thu nghiệm với hai thuật toán phía trên, các tâm Ξ được tạo bởi PDE Toolbox trong MATLAB

(phiên bản 2017b).

28

Bài toán 2.5.1 Cho phương trình Poisson ∆u = −3π2 sin πx sin πy sin πz trong miền Ω = [0, 1]3 với đều kiện biên Dirichlet u|∂Ω = 0. Nghiệm chính xác được đưa ra bởi u(x, y, z) = sin πx sin πy sin πz.

Sai số trung bình phương (RRMS) định nghĩa như sau

(u (ζ) − ˆu (ζ))2 (cid:80) ζ∈Ξint . (2.36) Ec = RRMS(u, ˆu, Ξint) := (cid:118) (cid:117) (cid:117) (cid:117) (cid:117) (cid:116) (u (ζ))2 (cid:80) ζ∈Ξint

Ngoài các lỗi được đưa ra trong bảng 2.1 và hình 2.1(a), Hình 2.1(b) so sánh mật độ của ma trận hệ thống. Mật độ của ma trận A ∈ Rn×n là số lượng trung bình của các mục nhập khác nhau trên mỗi hàng, được tính là nnz(A)/n, trong đó nnz(A) là tổng số mục nhập khác 0 trong A. Lưu ý rằng chúng tôi đo mật độ của ma trận (2.3) và của ma trận độ cứng phần tử hữu hạn sau khi loại bỏ điều kiện biên Dirichlet. Do

đó mật độ thấp hơn cho sự rời rạc thô trong đó tỉ lệ phần trăm đáng kể hơn của các nút nằm trên ranh giới. Mật độ đạt gần 16 đối với RBF-FD

1 và RBF-FD 2 đối với các nhóm nút lớn, trong khi đó đương nhiên là

RRMS lỗi trên các nút bên trong (Ec)

#Ξint

FEM

RBF-FD 1

RBF-FD 2

RBF-FD 3

33

9.19e-02

5.79e-02

4.87e-02

8.44e-02

80

7.19e-02

3.38e-02

1.80e-02

4.15e-02

179

4.64e-02

2.03e-02

1.25e-02

3.24e-02

479

2.19e-02

1.01e-02

7.21e-03

1.78e-02

1008

1.31e-02

4.47e-03

3.90e-03

1.34e-02

2213

7.58e-03

4.11e-03

2.34e-03

8.45e-03

4633

4.72e-03

2.49e-03

1.51e-03

4.83e-03

9684

3.12e-03

1.42e-03

9.93e-04

3.11e-03

19776

2.02e-03

8.16e-04

6.66e-04

1.85e-03

41409

1.06e-03

4.93e-04

3.42e-04

1.79e-03

Bảng 2.1: RRMS lỗi đối với nghiệm chính xác cho bài toán thử nghiệm số 2.5.1.

giống nhau đối với FEM và RBF-FD 3 và cách tiếp cận 14.

29

(b) Mật độ ma trận

(a) Lỗi trên các nút bên trong (Ec)

Hình 2.1: Bài toán kiểm tra 2.5.1. Biểu đồ đầu tiên trực quan hóa lỗi RRMS trên các nút bên

trong. Biểu đồ thứ hai cho thấy mật độ của ma trận hệ thống thưa thớt của hệ thống tuyến tính

(2.3) cho phương pháp RBF-FD và ma trận độ cứng cho FEM.

Bài toán 2.5.2 Phương trình Poisson ∆u = 3ex+y+z trong miền hình cầu đơn vị Ω = {(x, y, z) ∈ R3 : x2 + y2 + z2 ≤ 1} với các điều kiện biên Dirichlet được chọn sao cho nghiệm chính xác là u(x, y, z) = ex+y+z.

RRMS lỗi trên các nút bên trong (Ec)

#Ξint

FEM

RBF-FD 1

RBF-FD 2

RBF-FD 3

349

4.20e-03

1.40e-03

1.70e-03

1.77e-03

650

2.61e-03

1.11e-03

1.19e-03

8.68e-04

1318

1.66e-03

8.48e-04

7.74e-04

8.68e-04

2559

1.02e-03

4.90e-04

4.98e-04

4.93e-04

5254

6.04e-04

2.51e-04

2.18e-04

3.15e-04

10662

3.80e-04

1.72e-04

1.56e-04

1.80e-04

21777

2.38e-04

7.30e-05

8.46e-05

4.99e-04

43813

1.45e-04

4.67e-05

4.03e-05

1.43e-03

Bảng 2.2: RRMS lỗi đối với nghiệm chính xác cho bài toán thử nghiệm số 2.5.2.

Các kết quả được báo cáo trong Bảng 2.2 và hinh 2.2. Tam giác thô nhất thu được bằng cách tạo lưới với Hmax = 0.25 và có 349 đỉnh bên trong. Chúng tôi đánh giá hiệu suất của các thuật toán bằng cách sử dụng cùng một số đo Ec như trước đây.

30

(b) Mật độ ma trận

(a) Lỗi trên các nút bên trong (Ec)

Hình 2.2: Bài toán thử nghiệm 2.5.2. RRMS lỗi trên các nút bên trong và ma trận mật độ.

31

Kết luận

Trong quá trình tìm hiểu và thực hiện đề tài luận văn, trên cơ sở các kiến thức được học tập trong toàn khóa học và tìm hiểu tại các tài liệu

thuộc lĩnh vực của đề tài, em đã hoàn thành luận văn và đạt được các kết quả như sau:

- Trình bày những kiến thức cơ sở xung quanh luận văn. - Trình bày lược đồ phương pháp sai phân hữu hạn giải phương trình

đạo hàm riêng trong không gian ba chiều.

- Trình bày phương pháp nội suy hàm cơ sở bán kính, ứng dụng tính véctơ trọng số RBF.

- Trình bày cách xác định véctơ trọng số RBF trong không gian ba chiều. - Đưa ra hai bài toán thử nghiệm số để minh họa cho các Thuật toán

đã trình bày bên trên.

32

Tài liệu tham khảo

Tiếng Việt

[1] Tạ Văn Đĩnh (2002), “Phương pháp sai phân và phương pháp phần

tử hữu hạn”, Nhà xuất bản khoa học và kỹ thuật, Hà nội.

[2] Đặng Thị Oanh (2012), Phương pháp không lưới giải phương trình

Poisson, Luận án TS toán học, Viện công nghệ thông tin – Viện hàn lâm khoa học Việt Nam.

[3] Đàm Văn Mạnh (2013), “Nghiên cứu một số phương pháp nội suy

và xấp xỉ hàm số”, Luận văn thạc sỹ toán, Trường Đại học khoa học, Đại học Thái Nguyên, Thái Nguyên.

Tiếng Anh

[4] I. Boztosun, A. Charafi (2002), An analysis of the linear advection

– diffusion equation using mesh-free and mesh-dependent methods. Engineering Analysis with Boundary Elements, 26(1), 889–895.

[5] I. Boztosun, A. Charafi, D. Boztosun (2003), On the Numerical

Solution of Linear Advection-Diffusion Equation using Compactly Supported Radial Basis Functions, Meshfree Methods for Partial

Differential Equations, Lecture Notes in Computational Science and Engineering, 26, 63-73.

[6] M. D. Buhmann (2003), Radial Basis Functions, Cambridge Uni-

versity Press, New York, NY, USA.

33

[7] G. F. Fasshauer (2007), Meshfree Approximation Methods with

MATLAB, World Scientific Publishing Co., Inc., River Edge, NJ, USA.

[8] C. K. Lee, X. Liu, S.C. Fan (2003), Local multiquadric approxima- tion for solving boundary value problems. Comput. Mech. 30 (5–6),

396–409.

[9] A. I. Tolstykh and D. A. Shirobokov (2003), On using radial basis functions in a ‘finite difference mode with applications to elasticity

problems. Computational Mechanics, 33(1): 68-79.

[10] H. Wendland (2005), Scattered Data Approximation, Cambridge

University Press.

[11] G. B. Wright, Grady and B, Fornberg (2006), Scattered node com- pact finite difference-type formulas generated from radial basis func-

tions, J. Comput. Phys., 212(1): 99-123.

[12] Oleg Davydov and D. T. Oanh (2011), Adaptive meshless centres

and RBF stencils for Poisson equation. J. Comput, Phys., 230(1): 287-304.

[13] Oleg Davydov and D. T. Oanh (2011), On the optimal shape pa-

rameter for Gaussian radial basis function finite difference approx- imation of the Poisson equation, Computers and Mathematics with

Applications, 62(1): 2143-2161.

[14] Dang Thi Oanh, Oleg Davydov and Hoang Xuan Phu (2017), Adap-

tive RBF-FD method for elliptic problems with point singularities in

2D. Applied Mathematics and Computation, 313: 474-491.

[15] Oleg Davydov, Dang Thi Oanh and Ngo Manh Tuong (2019),

Octant-Based Stencil Selection for Meshless Finite Difference Meth- ods in 3D. Vietnam Journal of Mathematics.