Phương pháp số cho phương trình Helmholtz

Năm học 2016 - 2017 51 PHƯƠNG PHÁP SỐ CHO PHƯƠNG TRÌNH HELMHOLTZ Lê Thị Minh Thảo (Sinh viên năm 4, Khoa Toán – Tin học) GVHD: TS Nguyễn Thành Nhân 1. Mở đầu Được nghiên cứu cách đây chưa lâu, nhưng phương trình Helmholtz thu hút được sự nhiều sự quan tâm, chứng minh tính chính quy của nghiệm và nhiều cách giải được đưa ra nhằm tìm lời giải số cho phương trình này. Một trong những cách đó là dùng định lí Green để đưa bài toán về phương trình tích phân Lippmann – Schwinger. Trong nghiên cứu này, chúng ta không nhắc lại toàn bộ nhưng làm rõ hơn ở một số điểm về phương pháp số được đề cập đến trong bài báo của Vainikko, viết chương trình cho MATLAB nhằm tìm ra một số kết quả số minh họa, thông qua đó đưa ra một ý tưởng xấp xỉ khác và đánh giá sơ bộ về sai số, cho phương trình Lippmann – Schwinger, chính là dùng hệ đa thức Bernstein. 2. Nội dung Bài toán tán xạ cho phương trình Helmholtz trong môi trường không đồng tính: Giả sử miền không đồng tính là trơn hoặc trơn từng khúc; đồng thời có giá compact chứa gốc tọa độ. Chỉ số khúc xạ nb   , 2n  hoặc 3n  ,   1b x  bên ngoài miền không đồng tính; thỏa mãn b trơn hoặc trơn từng khúc. Bài toán được phát biểu như sau: Tìm : nu  ( 2n  hoặc 3n  ) thỏa mãn      2 0u x b x u x   , nx , (1) i su u u  , (2) 1 2lim 0 n s s r x ur i u r           đều, với  0,1x S x  (3) với 0  là số sóng. (3) được gọi là điều kiện bức xạ Sommerfeld. Có nhiều cách tiếp cận cho bài toán tán xạ này, ở đây ta nghiên cứu cách tiếp cận bằng đưa về phương trình tích phân. Áp dụng định lí Green, ta có hệ (1) – (3) suy ra phương trình tích phân Lippmann – Schwinger (tham khảo [2])          2 n iu x u x x y a y u y dy     (4) Kỉ yếu Hội nghị sinh viên NCKH 52 với 1a b  trơn hoặc trơn từng khúc, và có giá compact, với      10 , 24 1 , 3 4 ir i H r n r e n r         , 0r  , (5)  1 0H là hàm cầu Hankel loại một bậc 0 có công thức      1n n nH j z iy z  với  nj z là hàm cầu Bessel loại một và  ny z là hàm cầu Bessel loại hai được tìm ra theo       2 0 1 : 2 !1 3 2 2 1 p n p n p p t j t p n p         ,            2 1 0 2 ! 1 : 2 ! 2 ! 2 1 2 3 2 2 1 p p n n n p p n t y t n p n n n p               . Mệnh đề 0.1. Nếu  2 3u C là một nghiệm của phương trình Helmholtz với điều kiện Sommerfeld, khi đó u là nghiệm của (4). Ngược lại, nếu  3u C là nghiệm của (4) khi đó  2 3u C và u là một nghiệm của phương trình Helmholtz. Mệnh đề 0.2. Với mỗi 0k  phương trình Lippmann – Schwinger là giải được duy nhất nghiệm, hay nói cách khác, bài toán tán xạ là giải được duy nhất nghiệm. Sau đây, ta trình bày về nguyên lí chung của phương pháp số. Nguyên lí chung là, thay vì giải bài toán trong không gian vô hạn chiều, ta sẽ xấp xỉ không gian vô hạn chiều bởi một không gian hữu hạn chiều, nghiệm số chính là nghiệm có được khi giải bài toán trong không gian hữu hạn chiều. Quan trọng là không gian hữu hạn chiều sẽ rất gần với không gian vô hạn chiều khi số chiều dần ra vô cùng. Trong suốt báo cáo này, chúng ta sẽ thực hiện tìm nghiệm số theo phương pháp xấp xỉ sau: - Chọn một hệ cơ sở Hamel (hệ cơ sở đếm được) của hàm  2 RL G , RG là hình hộp tâm tại gốc tọa độ và có bán kính là R . - Xấp xỉ các hàm có mặt trong phương trình Lippman – Schwinger theo hệ cơ sở này dựa trên vector mốc nội suy. - Thiết lập hệ phương trình tuyến tính với ẩn là hệ số trong chuỗi hàm khai triển (dùng phương pháp đồng nhất). - Giải các hệ phương trình tuyến tính để tìm các hệ số trong chuỗi hàm khai triển. Năm học 2016 - 2017 53 2.1. Về xấp xỉ cho phương trình Lippman – Schwinger bằng chuỗi Fourier Đặt  : , 1, ,nR kG x x R k n     . Tham khảo phần 2 trong bài báo của Vainikko, dễ thấy thay vì xét (4) ta sẽ giải phương trình tích phân tuần hoàn đa chu kì, ẩn là      v x a x u x :             RG v x a x f x a x K x y v y dy   , với v au . (1.2) Với Nv là chuỗi Fourier hữu hạn của v , ta sẽ tìm nghiệm số xấp xỉ của (1.2) bằng cách giải phương trình xấp xỉ    N N N Nv Q af Q aKv  , (1.3) với K là toán tử tích phân:        RG Kv x K x y v y dy  có hạt nhân là hàm cut – off của  trên RG . Do đó, ma trận của phương pháp sắp xếp (1.3) được cho bởi NN N A v g , 1 NNN N N NA I a F K F   ,  N Ng af . (1.*) Trong đó, - NF và 1 NF  là các ma trận biến đổi Fourier thuận và ngược. - ,a f là các hàm số như đã cho trong phần mở đầu ở phương trình Helmholtz. - NK là ma trận đường chéo có đường chéo là giá trị của  K x tại các mốc nội suy là nút của lưới cơ sở. Định lí 1.4. Giả sử các hàm a và f thỏa mãn (6), và bài toán thuần nhất tương ứng với (1) – (3), với 1  , 0iu  , có duy nhất nghiệm tầm thường. Khi đó phương trình (1.2) có nghiệm duy nhất v H  , phương trình sắp xếp (1.3) có nghiệm duy nhất Nv với 0N N , và N Nv v c v Q v c v N         , 0    . (1.9) 2.2. Mở rộng cho phương pháp sắp xếp (collocation method) Cũng nhằm tìm ra nghiệm số của phương trình Lippmann – Schwinger và dùng ý tưởng xấp xỉ, ta chọn hệ cơ sở là các hàm trong tập 2: exp , , n h h x jhX mh j Dh            , Kỉ yếu Hội nghị sinh viên NCKH 54 với h là miền chứa giá của  q x , h  , tham số D là số thực dương cố định và h là tham số rời rạc. Ta có một số kết quả sau (tham khảo các bổ đề và chứng minh trong [4]) Phương trình tích phân (2.5) được giải bằng phương pháp sắp xếp: Tìm h hu X sao cho             n h h ku mh q mh Ku mh q mh mh y g y dy     (2.7) với mọi điểm trên lưới hmh . Do đó, các hệ số  mu của nghiệm rời rạc     2 2/ h x mh Dh h m mh u x u e     được tính từ hệ phương trình tuyến tính         nh m m j j k jh u q mh a u q mh mh y g y dy        (2.8) với hmh , trong đó  ja jh với     2 2/ n y Dh kx x y e dy    (2.9) Gọi hQ là phép chiếu nội suy lên không gian các Gaussian đã định nghĩa phía trước. Như vậy ta chỉ còn đánh giá hạng tử thứ hai  hK I Q q . Bổ đề 2.5. Giả sử hàm liên tục u thỏa mãn    : 1 n N N u Fu t t dt          , với  thỏa mãn 0 1  . Khi đó sai số nội suy có thể được đánh giá như sau           0 0 2 2 ! N N h N u xDhu x Q u x c h u a x                   , với hằng số c nào đó không phụ thuộc u . Định lí 2.6. Giả sử nghiệm  của phương trình (2.1) thỏa mãn tính trơn của điều kiện của kết quả phía trên với 2 2N m  và công thức thể tích cho vế phải của (2.5) được tạo ra bởi 2m hoặc 2m . Khi đó sai số giữa h so với  của (2.8) thỏa mãn Năm học 2016 - 2017 55       21 N h ux x c Dh c h     . Kế tiếp ta đưa ra mở rộng bằng cách dùng chuỗi đa thức Bernstein để xấp xỉ. 2.3. Dùng hệ đa thức Bernstein để xấp xỉ phương trình Lippmann - Schwinger Định nghĩa 3.1. Cho 0 n  , 0  và các số  1, , 1, ,mk k n  là các chỉ số khi 0  . Khi đó ta có định nghĩa đa thức Bernstein cấp  như sau      00 0 1 1 2 1 1 2 : , , , , , , , , , 1 m m m k n k k k B f x f x x x x x                                    với  1: , , mk k     ,  0 1 2: 0, , ,0, , , , , ,0m        , và thỏa mãn j nằm ở vị trí thứ jk . Khi 2n  ta có dạng tường minh của đa thức Bernstein như sau:       1 2 1 1 2 21 2 1 2 1 2 1 21 2 , , 1 2 1 1 2 2 0 0 1 21 2 , : , 1 1 m m m k m kk k f m m k k m mk kB x x f x x x x k km m                    . Ta cũng có một cách viết khác tường minh hơn cho đa thức Bernstein:       1 1 , , , 1 0 11 1, , , , : , , 1 j jj n j j n m kj kn f m m n j j k m jjn j n mkkB x x f x x km m                           . Định lí 3.4. Cho  : : 0,1 nf I   liên tục thỏa điều kiện Lipschitz     22f x f y L x y   trên I , khi đó ta có bất đẳng thức     1 1 2 , , , 2 1 1 2n n f m m j j LB x f x m            . Ta sẽ giải phương trình Lippmann – Schwinger bằng phương pháp xấp xỉ đa thức Bernstein    N N N NB af B aK   , trong đó K là toán tử tích phân        RG Ku x x y u y dy   . Kỉ yếu Hội nghị sinh viên NCKH 56 Thực hiện tương tự như phần 3.4 của [1], với lưu ý là NB là ma trận biến đổi các hệ số đa thức Bernstein (theo công thức của đa thức Bernstein, chính là ma trận đơn vị), khi đó dạng ma trận của phương pháp này được cho bởi hệ phương trình tuyến tính N N NA g  , 1 N NNN N NA I a K  B B ,  NN Ng af B . (3.5) Lưu ý NK là ma trận đường chéo trong đó đường chéo là vector có tọa độ là các hệ số Bernstein của hàm  tại các mốc nội suy là điểm trên lưới. Định lí 3.6. Giả sử ,a f là các hàm số trên không gian n thỏa mãn  supp 0,a B  ,  ,2 na W  ,  ,2loc nf W  với 2 n   , và bài toán tán xạ thuần nhất với 1  , 0iu  chỉ có duy nhất nghiệm tầm thường. Khi đó phương trình             RG x a x f x a x x y y dy     (3.6) có duy nhất nghiệm trong  2 nC , và phương trình xấp xỉ có nghiệm duy nhất  nN N x  và ta có bất đẳng thức 1 2 1 2 22 2N N c B c N        . (3.7) 2.4. Tốc độ hội tụ Định lí 4.4. Cho D X , với X là không gian Banach và :f X Y , Y là không gian Banach,  pf C D . Nếu  pf có module liên tục  , khi đó tổng riêng của chuỗi Fourier của f hội tụ về f với tốc độ      ln 2N p Nf x S f x K N N          . (4.1) Nếu f thỏa mãn điều kiện  - Holder, khi đó      lnN Nf x S f x K N    . (4.2) Mệnh đề 4.5. Cho f liên tục trên  0,1 và ,n fB là đa thức Bernstein ứng với f . Theo kết quả quen thuộc thì ,n fB hội tụ đều về f trên  0,1 và tốc độ hội tụ của ,n fB thỏa bất đẳng thức  ,n f fB f c n  , (4.3) Năm học 2016 - 2017 57 với 5 4 c  và f là module liên tục của f trên  0,1 . Nếu giá trị của c trong (4.3) được thay bởi 4306 837 6 5832  thì (4.3) trở thành đánh giá tốt nhất, nghĩa là không tồn tại giá trị c hay lũy thừa nào của n làm cho bất đẳng thức (4.3) chặt hơn. Dù các đánh giá về tốc độ hội tụ phía trên là cho các hàm một biến, ta cũng có được cái nhìn tổng quan rằng tốc độ hội tụ của chuỗi Fourier “nhanh” hơn so với đa thức Bernstein rất nhiều, nhất là khi  càng lớn. Do đó dù ma trận biến đổi hệ số của đa thức Bernstein đơn giản hơn so với ma trận biến đổi Fourier , người ta vẫn hay dùng chuỗi Fourier để xấp xỉ. 2.5. Chương trình MATLAB thực hiện phương pháp sắp xếp (collocation method) Chương trình được thực hiện qua 9 bước 1. Chuyển hàm a và f cho trước thành các function tương ứng trong MATLAB. 2. Viết function xây dựng ma trận biểu diễn các nút của một lưới có độ mịn h trên hình vuông tâm tại gốc tọa độ, cạnh 2R , (để tiện lợi ta gọi là hình vuông cơ sở), với 2Rh N  . 3. Viết function tính vector có tọa độ là các giá trị của một hàm tại các điểm là nút của lưới cơ sở. 4. Viết function tính ma trận biến đổi Fourier hai chiều. 5. Viết function tính vector có tọa độ là các giá trị của hàm cut – off  K x tại các điểm là nút của lưới cơ sở. 6. Viết function chuyển đổi vector thành ma trận đường chéo mà đường chéo này chính là tọa độ trong vector đã cho (theo thứ tự). 7. Viết function tính vector Nv có tọa độ là giá trị của hàm cần tìm tại các nút của lưới. 8. Viết function tính vector giá trị của chuỗi Fourier của Nv tại tất cả các điểm trên hình vuông cơ sở. 9. Biểu diễn hàm số tìm được dưới dạng đồ thị ba chiều. Do ta chỉ quan tâm đến phần thực của vector Nv nên hình vẽ dưới đây biểu diễn phần thực của chuỗi khai triển Fourier của hàm cần tìm. Kỉ yếu Hội nghị sinh viên NCKH 58 Hình 1. 3, 3N R  . Hình 2. 7, 3N R  Năm học 2016 - 2017 59 Hình 3. 15, 4N R  . Hình 4. 31, 4N R  . Kỉ yếu Hội nghị sinh viên NCKH 60 Các hình vẽ trên cho ta cái nhìn trực quan về hàm cần tìm – chính là hàm sóng. Khi cho N càng lớn thì đồ thị càng cho ta hình dạng rõ ràng của sóng hơn, càng ra xa biên độ càng giảm (do sự tán xạ) và độ chính xác càng cao (thể hiện rõ sự tán xạ). 3. Kết luận Bài báo cáo này đã trình bày một phương pháp tìm nghiệm số xấp xỉ cho phương trình Helmholtz thông qua phương trình tích phân Lippmann – Schwinger. Cả ba phương pháp đưa ra đề dựa trên ý tưởng là xấp xỉ một hàm khả tích (hoặc trên không gian Sobolev) bằng những hàm quen thuộc (chuỗi Fourier, đa thức) và được trình bày những đánh giá sai số. Dù trên lí thuyết là có thể tìm nghiệm số theo những cách này tuy nhiên ta nhận thấy rằng cần cải tiến thuật toán để đảm bảo thời gian chạy thuật toán. Trong báo cáo này, ngoài việc chuyển ngữ và biên tập lại một số kết quả trong [1] và [4], chúng tôi đã thực hiện một số nội dung bổ sung: 1. Bổ sung một số chứng minh cho các mệnh đề trong [1]; 2. Làm rõ hai thuật toán DFT và FFT sửa dụng trong [1], đồng thời giải thích một số cách cải tiến dựa trên DFT và FFT; 3. Triển khai và chứng minh ý tưởng dùng đa thức Bernstein để tìm nghiệm số xấp xỉ của bài toán (toàn bộ chương 3); 4. Trình bày sơ lược về tốc độ hội tụ để đưa ra nhận xét nên dùng hệ cơ sở nào khi xấp xỉ. Tuy nhiên, do năng lực và thời gian có hạn, trừ đoạn chương trình cho MatLab ở chương 1 dành cho thuật toán của Vainikko, hai phần sau vẫn chưa có đoạn chương trình. Đồng thời, chúng tôi cũng chưa đánh giá được thời gian tiêu tốn cho hai thuật toán ở chương 2 (ở chương 1, trong [1] đã có đưa ra tính toán). Quan trọng hơn, ở chương 2 và 3 chúng tôi chỉ dừng lại ở việc chứng minh tính đúng đắn nhưng vẫn chưa đưa ra một thuật toán và ví dụ tường minh. Chúng tôi vẫn đang trong quá trình hoàn thiện. Về tương lai, phương trình Helmholtz vẫn còn đang và sẽ được nghiên cứu mạnh mẽ, phương pháp xấp xỉ như trên vẫn đang đòi hỏi sự cải tiến và hơn nữa phương trình này vẫn đang đòi hỏi những ý tưởng và lời giải mới. Năm học 2016 - 2017 61 TÀI LIỆU THAM KHẢO 1. Gennadi Vainikko (2000), Fast solvers of the Lippmann – Schwinger equation, In P. P. Gilbert, J. Kajiwara, and Y. S. Xu, Editors, Direct and Inverse Problems of Mathematical Physics (Newark, DE, 1997), pages 423 – 440. Kluwer, Dordrecht. 2. Stephen Roberts, Lecture 7 – The Discrete Fourier Transform, University of Oxford. 3. D. Colton and R. Kress (1992), Inverse Acoustic and Electromagnetic Scattering Theory, Springer. 4. F. Lanzara, V. Maz’ya, and G. Schmidt, Numerical Solution of the Lippmann – Schwinger Equation by Approximate Approximations, The Journal of Fourier Analysis and Applications, “Online First”. 5. Maz’ya, V. and Schmidt, G. (1995), Approximate approximations and the cubature of potentials, Rend. Mat. Acc. Lincei, 6(9), 161 – 164. 6. Powell, M. J. D (1992), The theory of radial basis functions in 1990, in Advance in Numerical Analysis, Vol 2: Wavelets, Subdivision Algorithm, and Radical Basis Functions, Light, W. Ed, 105 – 210, Clarendon Press, Oxford. 7. Walter Rudin, Principle of Mathematical Analysis, Third Edition, McGraw Hill Inc. 8. Adrian Fellhauer (2016), Approximation of smooth functions using Berstein polynomials in multiple variables, arXiv: 1609.01940v1. 9. Clemens Hitzinger (1974), Dissertation Simulation and Inverse Modeling of Semiconductor Manufacturing Process: 7.4 Multivariate Bernstein Polynomials, Luftbadgasse 11, A-1060 Wien, Matrikelnummer e9425899, geboren am 9. in Linz. 10. Jihad Titi, Jurgen Garloff (2010), Matrix methods for the Bernstein Form and their application in Global Optimization, Presentation. 11. Dunham Jackson (1930), The theory of Approximation, American Mathematical Scociety Colloquium Publications, Volume 11. 12. P. C. Sikkema, Der Wert (1961), Konstanten in der Thoerie der Approximation mit Bernstein – Polynomen, Numer. Math, 107 – 116.