Phương trình tích phân

3 MỞ ĐẦU Phương trình tích phân là một trong những loại phương trình quan trọng, xuất phát từ nhiều bài tốn ứng dụng. Đặc biệt đối với nhiều phương trình đạo hàm riêng mà bằng cách dùng hàm Green, người ta quy về việc giải một phương trình tích phân dạng ( ) ( ) ( ) ( ) D f t K t s f s ds g tα + − =∫ , (1) hay ( ) ( ) ( ) D K t s f s ds g t− =∫ , (2) với t D′∈ , trong đĩ nhân K và hàm g cho trước và f là Nn hàm cần tìm. Tổng quát, (1) được gọi là phương trình tích phân Fredholm loại hai và là bài tốn chỉnh (với các mêtríc thơng thường trên các khơng gian hàm) và (2) là phương trình loại một và là bài tốn khơng chỉnh. Trong trường hợp D D′= = ℝ hay nℝ , người ta cĩ thể dùng phép biến đổi Fourier để giải (1) cũng như xây dựng nghiệm chỉnh hĩa tường minh cho (2). Nội dung luận văn nhằm tìm một cách tiếp cận khác cho các trường hợp D ≠ ℝ hay nℝ , trong đĩ các phương trình (1) và (2) được xấp xỉ trực tiếp bằng các hệ phương trình tuyến tính với Nn cần tìm chính là giá trị của Nn hàm tại các điểm nút trong miền xác định của nĩ. Hơn nữa, với một số các điều kiện thích hợp, hệ phương trình tuyến tính nhận được cĩ ma trận các hệ số thuộc loại đặc biệt : ma trận Toeplitz, ma trận Toeplitz đối xứng hay ma trận Circulant Toeplitz mà để giải nĩ, luận văn xây dựng một giải thuật đệ quy nhằm xác định nghiệm. Giải thuật này cĩ ưu điểm về mặt tính tốn do nĩ cĩ độ phức tạp nhỏ hơn hẳn các giải thuật giải hệ phương trình tuyến tính thơng thường. 4 Với nội dung chính nêu trên, luận văn bao gồm 3 chương, trong đĩ Chương 1 giới thiệu một số khái niệm và kết quả về các loại ma trận Toeplitz, sự liên hệ giữa ma trận Toeplitz với phương trình tích phân dạng (1) và (2). Đặc biệt, một số kết quả về giá trị riêng cũng như vectơ riêng của các loại ma trận này cũng được khảo sát. Điều này sẽ giúp ích sau này trong việc kiểm sốt các sai số tính tốn khi giải số bằng máy tính. Chương 2 nhằm xây dựng một số kết quả mà đặc biệt là xây dựng giải thuật cũng như chương trình giải xấp xỉ phương trình tích phân Fredholm loại hai. Chương 3 đưa ra giải thuật đệ quy giải hệ phương trình tuyến tính cĩ ma trận các hệ số là ma trận loại Toeplitz. Từ đĩ, tính tốn cũng như so sánh trên một số ví dụ cụ thể cho các hệ phương trình tuyến tính cũng như một số phương trình tích phân cụ thể. 5 Chương 1 MA TRẬN TOEPLITZ & CIRCULANT TOEPLITZ Trong chương này, chúng tơi trình bày một số kết quả cho ma trận Toeplitz, ma trận Toeplitz đối xứng và ma trận Circulant Toeplitz, trong đĩ ma trận loại này xuất hiện khi ta xấp xỉ một phương trình tích phân loại Fredholm với nhân dạng tích chập về một hệ phương trình tuyến tính. 1. MA TRẬN TOEPLITZ 1.1. Khái niệm ma trận Toeplitz Một ma trận vuơng cấp n được gọi là ma trận Toeplitz (hay ma trận hằng đường chéo), ký hiệu n k, j k, j 0,1, ,n 1T t = − =   … , khi k, j k jt t −= , với mọi k, j 0,1, ,n 1= −… , nghĩa là ma trận cĩ dạng 0 1 2 (n 1) 1 0 1 (n 2) n 2 1 0 (n 3) n 1 n 2 n 3 0 t t t t t t t t T t t t t t t t t − − − − − − − − − − − −         =          … … … ⋮ ⋮ ⋮ ⋱ ⋮ … (1.1) Loại ma trận này xuất hiện nhiều trong các bài tốn ứng dụng. Chẳng hạn, xét phương trình Fredholm loại một, ( ) ( ) ( ) b a g t K t,s f s ds= ∫ (1.2) với ( )f t là Nn hàm cần tìm, ( )g t và ( )K t,s là các hàm cho trước, trong đĩ ( )K t,s được gọi là nhân của phương trình tích phân. 6 Giả sử a, b hữu hạn. Ta chia đoạn [ ]a ;b thành n 1− đoạn bằng nhau [ ]0 1s ;s , [ ]1 2s ;s , ..., [ ]n 2 n 1s ;s− − , mỗi đoạn cĩ độ dài b ah n 1 − = − , bằng các điểm chia is a ih= + , i 0, n 1= − . Ta cĩ thể xấp xỉ tích phân ở vế phải bằng cơng thức cầu phương, ( ) ( ) b n 1 i i i 0a y s ds w y s − = =∑∫ trong đĩ ( )iw là những trọng số của cơng thức cầu phương. Trường hợp đơn giản nhất, ta chọn iw 1= . Áp dụng cơng thức trên vào phương trình (1.2), ta được ( ) ( ) ( )n 1 i i i 0 g t K t,s f s , t [a ;b] − = = ∀ ∈∑ (1.3) Khi đĩ, tại các điểm jt t= , với j 0, n 1= − , ta cĩ hệ phương trình ( ) ( ) ( )n 1j j i i i 0 g t K t ,s f s − = =∑ và tại những điểm j js t= , ta được hệ phương trình tuyến tính theo giá trị của Nn hàm f tại các điểm is , ( ) ( ) ( )n 1j j i i i 0 g s K s ,s f s − = =∑ (1.4) Trong trường hợp nhân K cĩ dạng ( ) ( )K t,s t s= ϕ − , trong đĩ ϕ là một hàm số thực cho trước, ta cĩ ( )j, i j i j it s s t −= ϕ − = , và với ( )j jg s g= , ( )i if s f= , phương trình (1.4) được viết lại thành n 1 n 1 j j, i i j i i i 0 i 0 g t f t f − − − = = = =∑ ∑ 7 nghĩa là ta được hệ o 0 1 1 2 2 (n 1) n 1 0 1 0 0 1 1 2 (n 2) n 1 1 2 0 1 1 0 2 (n 3) n 1 2 n 1 0 n 2 1 n 3 2 0 n 1 n 1 t f t f t f t f g t f t f t f t f g t f t f t f t f g t f t f t f t f g − − − − − − − − − − − − − − − − − + + + + =  + + + + =  + + + + =    + + + + = … … … …… … (1.5) Hệ phương trình (1.5) được viết lại dưới dạng phương trình ma trận nT f g= , trong đĩ 0 1 2 (n 1) 1 0 1 (n 2) n 2 1 0 (n 3) n 1 n 2 n 3 0 t t t t t t t t T t t t t t t t t − − − − − − − − − − − −         =          … … … ⋮ ⋮ ⋮ ⋱ ⋮ … là ma trận Toeplitz cấp n và 0 1 n 1 g g g g −       =       ⋮ là những ma trận đã biết, và 0 1 n 1 f f f f −       =       ⋮ là ma trận Nn cần tìm. 8 Một trường hợp đặc biệt của ma trận Toeplitz là ma trận Toeplitz đối xứng dạng 0 1 2 n 1 1 0 1 n 2 n 2 1 0 n 3 n 1 n 2 n 3 0 t t t t t t t t T t t t t t t t t − − − − − −        =         … … … ⋮ ⋮ ⋮ ⋱ ⋮ … . Chú ý rằng ma trận Toeplitz đối xứng được hồn tồn xác định chỉ bằng n số thực 0 1 n 1t , t , t − ∈ℝ . 1.2. Vectơ riêng và giá trị riêng Nhắc lại rằng số phức α được gọi là một trị riêng của ma trận A nếu tồn tại vectơ x khác vectơ khơng sao cho Ax x= α (1.6) và khi đĩ vectơ x được gọi là một vectơ riêng tương ứng với trị riêng α . Trường hợp đặc biệt khi A là ma trận Hermite, nghĩa là *A A= , trong đĩ *A là ma trận liên hợp chuyển vị của A thì các trị riêng của A là số thực. Khi đĩ, ta cĩ thể giả sử dãy giá trị riêng ( )iα của ma trận A là dãy khơng tăng, nghĩa là 0 1 2α ≥ α ≥ α ≥… Xét dãy ma trận Toeplitz Hermite n k j k, j 0,1,...,n 1T t − = − =   ứng với tập giá trị riêng { }n,i ; i 0,1,2,..., n 1τ = − . Dáng điệu tiệm cận của các giá trị riêng { }n,i ; i 0,1,2,...,n 1τ = − này được cho bởi Định lý Szegư với điều kiện tồn tại sự liên hệ của chuỗi Fourier với các hệ số kt , ( ) [ ]ikk k f t e , 0;2 +∞ λ =−∞ λ = λ ∈ pi∑ (1.7) và 9 ( ) 2 ik k 0 1 t f e d 2 pi − λ = λ λ pi ∫ . (1.8) Điều kiện này cĩ nghĩa rằng dãy ( )kt xác định một hàm f và ngược lại, các kt được xác định từ một hàm số f. Khi đĩ, dãy các ma trận Toeplitz được ký hiệu bởi ( )nT f . Khi ( )nT f là ma trận Hermite thì *k kt t− = và f là hàm nhận giá trị thực. Khi đĩ với mỗi hàm F liên tục trên miền xác định của f , định lý Szegư cho ( ) ( )( )2n 1 n,k n k 0 0 1 1lim F F f d n 2 pi − →∞ = τ = λ λ pi ∑ ∫ (1.9) Chẳng hạn, với ( )F x x= , ta cĩ ( ) 2n 1 n,k n k 0 0 1 1lim f d n 2 pi − →∞ = τ = λ λ pi ∑ ∫ (1.10) và điều này cĩ nghĩa là trung bình cộng các trị riêng của nT (f ) hội tụ về tích phân của hàm f . Mặt khác, do vết ( )Tr A là tổng các phần tử trên đường chéo của A và cũng chính là tổng của các giá trị riêng của A khi A là ma trận Hermite, nên từ (1.10), ta suy ra ( )( ) ( )2n n 0 1 1lim Tr T f f d n 2 pi →∞ = λ λ pi ∫ (1.11) Tương tự, với sF(x) x= , ta được ( ) 2n 1 s s n,k n k 0 0 1 1lim f d n 2 pi − →∞ = τ = λ λ pi ∑ ∫ (1.12) Nếu f là hàm thực và các giá trị riêng thỏa n,k m 0 ; n,kτ ≥ > ∀ thì do hàm F(x) ln x= liên tục trên [ )m ; + ∞ , ta suy ra 10 ( ) 2n 1 n,i n i 0 0 1 1lim ln ln f d n 2 pi − →∞ = τ = λ λ pi ∑ ∫ (1.13) Mặt khác, do định thức của các ( )nT f chính là tích các giá trị riêng của nĩ, ( )( ) n 1n n,i i 0 det T f − = = τ∏ , (1.13) trở thành ( )( )( ) ( )21 n 1nn n,i n x i 0 0 1 1lim ln det T f lim ln ln f d n 2 pi − →∞ →∞ =   = τ = λ λ  pi  ∑ ∫ (1.14) Nếu f bị chặn dưới thì tất cả các giá trị riêng của nT (f ) sẽ bị chặn dưới và (1.14) mơ tả giới hạn của dãy định thức Toeplitz. 2. MA TRẬN CIRCULANT TOEPLITZ 2.1. Khái niệm ma trận Circulant Toeplitz Ma trận C được gọi là ma trận Cicurlant Toeplitz nếu nĩ là ma trận Toeplitz cĩ dạng 0 1 2 n 1 n 1 0 1 2 n 1 0 1 1 n 1 0 c c c c c c c c C c c c c c c − − − −        =         … ⋮ ⋱ ⋮ ⋱ ⋱ … (1.15) ký hiệu ( )k, j j k mod nC C c − = =   . Đây là một trường hợp đặc biệt của ma trận Toeplitz, trong đĩ là mỗi hàng dưới của ma trận là một dịch chuyển tuần hồn của hàng trên. Ma trận Cicurlant Toeplitz thường dùng để xấp xỉ cũng như biểu diễn dáng điệu tiệm cận của các ma trận Toeplitz. 11 2.2. Vectơ riêng và giá trị riêng Giá trị riêng kψ và vectơ riêng (k )y của C là nghiệm phương trình Cy y= ψ (1.16) hoặc tương đương với n phương trình sai phân m 1 n 1 n m k k k m k m k 0 k m c y c y y ; m 0, 1, ,n 1 − − − + − = = + = ψ = −∑ ∑ … . (1.17) Đổi chỉ số, ta được n 1 m n 1 k k m k k (n m) m k 0 k n m c y c y y ; m 0, 1, ,n 1 − − − + − − = = − + = ψ = −∑ ∑ … (1.18) Ta cĩ thể giải các phương trình sai phân này với cùng một phương pháp giải các phương trình vi phân bằng cách tìm nghiệm đặc biệt của nĩ. Do phương trình là tuyến tính với các hệ số hằng, một nghiệm đặc biệt hợp lý là kky = ρ (tương tự như dạng nghiệm mũ ( ) sty t e= trong phương trình vi phân tuyến tính). Thay vào (1.18), ta được n 1 m n 1 k m k (n m) m k k k 0 k n m c c − − − + − − = = − ρ + ρ = ψρ∑ ∑ hay n 1 m n 1 k n k k k k 0 k n m c c − − − − = = − ρ + ρ ρ = ψ∑ ∑ . Nếu ta chọn n 1−ρ = , nghĩa là ρ là một căn bậc n phức của đơn vị, ta nhận được giá trị riêng n 1 k k k 0 c − = ψ = ρ∑ (1.19) với vectơ riêng tương ứng 12 ( )T2 n 11y 1, , , , n − = ρ ρ ρ… (1.20) Cụ thể, với 2i m n m e pi −ρ = , ta cĩ giá trị riêng 2 imkn 1 n m k k 0 c e pi − − = ψ =∑ (1.21) và vectơ riêng tương ứng ( )( )2 i n 12 imn n T(m) 1y 1,e , ,e n pi −pi − − = … Từ định nghĩa của trị riêng và vectơ riêng, (m) (m) mCy y ; m 0,1, ,n 1= ψ = −… (1.22) và do (1.21) chính là biểu thức Fourier rời rạc (DFT) của dãy k(c ) , ta cĩ thể phục hồi k(c ) từ các kψ bằng cơng thức biến đổi Fourier ngược. Đặc biệt, ( ) ( ) 2 imk n 2 i l k m n n 1 n 1 n 1 2 ilm 2 ilm m k m 0 m 0 k 0 n 1 n 1 k l k 0 m 0 1 1 e c e e n n 1 c e c n pi pi − − − − −pi pi = = = − − = = ψ = = = ∑ ∑∑ ∑ ∑ (1.23) trong đĩ ta dùng tính chất 2 imk n n l k mod n m 0 n khi k mod n 0 e n 0 khi k mod n 0 pi − = = = δ =  ≠ ∑ (1.24) với ký hiệu Kronecker m 1 khi m 0 0 khi m 0 =δ =  ≠ Như vậy, giá trị riêng của ma trận Circulant là biến đổi Fourier rời rạc của dịng đầu tiên của nĩ, và ngược lại, dịng đầu tiên của ma trận Circulant là biến đổi Fourier rời rạc ngược của các giá trị riêng của nĩ. 13 Phương trình (1.22) cĩ thể được viết lại dưới dạng ma trận CU U= Ψ (1.25) trong đĩ (0) (1) (2) (n 1) 2 imk n m,k 0,1, , n 1 U y y y y 1 e n − pi − = −  =     =     … … là ma trận mà mỗi vectơ cột là một vectơ riêng và ( )kdiagΨ = ψ là ma trận chéo với các phần tử chéo là 0ψ , 1ψ , ..., n 1−ψ . Hơn nữa (1.24) cho thấy U là ma trận unita. Thật vậy, ký hiệu phần tử hàng k cột j của *UU bởi k, ja và chú ý rằng k, ja chính là tích của hàng thứ k của U, 2 imk n 1 e ;m 0,1, ,n 1 n pi −  = −    … với cột thứ j của *U , 2 imj n 1 e ;m 0,1, ,n 1 n pi  = −    … . Do đĩ ( ) ( ) 2 im j kn 1 n k, j k j mod n m 0 1 a e n pi − − − = = = δ∑ và ta cĩ *UU I= . Tương tự *U U I= . Do vậy (1.25) trở thành *C U U= Ψ (1.26) *U CUΨ = (1.27) 2.3. Các phép tốn trên ma trận Circulant Từ các kết quả nhận được trong phần trên, ta tổng kết một số tính chất quan trọng về phép tốn trên ma trận Circulant như sau : 14 Định lý. Mỗi ma trận Circulant C cĩ vectơ riêng ( )( )2 i n 12 imn n T(m) 1y 1,e , ,e n pi −pi − − = … , với m 0, 1, ,n 1= −… , tương ứng với các giá trị riêng 2 imkn 1 n m k k 0 c e pi − − = ψ =∑ và cĩ thể biểu thị dưới dạng *C U U= Ψ trong đĩ U cĩ các vectơ cột chính là các vectơ riêng và Ψ là ma trận chéo ( )kdiagΨ = ψ . Đặc biệt, mọi ma trận Circulant cĩ cùng các vectơ riêng cĩ thể biểu diễn bởi ma trận U và bất kỳ ma trận C cĩ dạng *C U U= Ψ là ma trận Circulant. Cho ( )k jC c −= và ( )k jB b −= là hai ma trận Circulant cấp n với các giá trị riêng tương ứng 2 imkn 1 n m k k 0 c e pi − − = ψ =∑ , 2 imkn 1 n m k k 0 b e pi − − = β =∑ . Khi đĩ (1) C và B giao hốn, *CB BC U U= = γ , trong đĩ ( )m mdiagγ = ψ β , và CB cũng là ma trận Circulant. (2) C B+ là ma trận Circulant và *C B U U+ = Ω trong đĩ ( )( )m m k m−Ω = ψ + β δ . (3) Nếu m 0 ; m 0, 1, ,n 1ψ ≠ = −… thì C là khơng suy biến và 1 1 *C U U− −= Ψ . 15 Chứng minh. Ta cĩ *C U U= Ψ và *B U U= Φ , với ( )mdiagΨ = ψ và ( )mdiagΦ = β . (1) Ta cĩ * * * *CB U U U U U U U U BC= Ψ Φ = ΨΦ = ΦΨ = và vì ΦΨ là ma trận chéo nên ta suy ra CB là ma trận Circulant. (2) *C B U( )U+ = Ψ + Φ (3) Nếu Ψ khơng suy biến, thì 1 * * 1 * 1 * *CU U U U U U U U UU I− − −Ψ = Ψ Ψ = ΨΨ = = . Do đĩ, 1 1 *C U U− −= Ψ . 16 Chương 2 PHƯƠNG TRÌNH FREDHOLM LOẠI HAI Trong chương này, chúng ta khảo sát phương trình Fredholm loại hai dạng ( ) ( ) ( ) ( ) b a f t K t,s f s ds g t= λ +∫ (2.1) với ( )f t là Nn hàm cần tìm, ( )g t và ( )K t,s là các hàm cho trước, trong đĩ ( )K t,s được gọi là nhân của phương trình tích phân. Khi ( )g t 0= phương trình tương ứng được gọi là thuần nhất. Ta khảo sát 3 bài tốn sau : Bài tốn 1 : Tìm nghiệm của phương trình khơng thuần nhất (2.1) với λ và ( )g t cho trước. Bài tốn 2 : Tìm giá trị riêng và vectơ riêng của nhân ( )K t,s , tức là tìm tham số λ , để phương trình thuần nhất ( ) ( ) ( ) b a f t K t,s f s ds= λ∫ cĩ nghiệm khơng tầm thường ( )f t . Giá trị λ nhận được gọi là giá trị riêng của ( )K t,s và hàm ( )f t 0≠ tương ứng được gọi là hàm riêng (vectơ riêng) của ( )K t,s . Bài tốn 3 : Các giải thuật giải số phương trình Fredholm loại hai và thực hành một số ví dụ cụ thể. 17 1. TÌM NGHIỆM XẤP XỈ CỦA PHƯƠNG TRÌNH FREDHOLM LOẠI HAI KHƠNG THUẦN NHẤT 1.1. Phương pháp Nystrom (thay nhân bằng tổng hữu hạn) Dùng cơng thức cầu phương ( ) ( )b n j j j 1a y s ds w y s = =∑∫ (2.2) trong đĩ ( )jw là những trọng số của cơng thức cầu phương và [ ]1 2 ns ,s ,...,s a ; b∈ . Thay cơng thức (2.2) vào (2.1), ta được phương trình xấp xỉ ( ) ( ) ( ) ( )n j j j j 1 f t w K t,s f s g t = = λ +∑ (2.3) Khi đĩ, với it t= , i 1,..., n= , ta cĩ phương trình ( ) ( ) ( ) ( )ni j i j j i j 1 f t w K t ,s f s g t = = λ +∑ ; i 1,..., n= (2.4) Đặt ( )i if f t= , ( )i ig g t= , ( )ij i jK K t ,s= và ij i j jK K w=ɶ thì (2.4) được viết lại thành hệ phương trình tuyến tính 1 11 1 12 2 1n n 1 2 21 1 22 2 2n n 2 n n1 1 n2 2 nn n n f K f K f ... K f g f K f K f ... K f g ... f K f K f ... K f g  = λ + λ + + λ +  = λ + λ + + λ +    = λ + λ + + λ + ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ (2.5) hay dưới dạng phương trình ma trận ( )1 K f g− λ ⋅ =ɶ . (2.6) Hệ (2.5) là hệ gồm n phương trình tuyến tính, n Nn. Giải hệ phương trình này, ta được giá trị của hàm f tại những điểm ( )it . Tại những điểm t khác, ta dùng 18 (2.3) làm cơng thức nội suy. Như vậy, biểu thức giải tích của nghiệm gần đúng cĩ thể lấy ở cơng thức (2.3). Để tìm các giá trị riêng gần đúng, trong (2.6) cho ig 0= , ta được hệ phương trình thuần nhất. Để hệ cĩ nghiệm khơng tầm thường thì định thức của hệ phải bằng khơng. Từ đĩ, ta được phương trình bậc n của λ . Giải phương trình này, ta được các giá trị riêng gần đúng 1 n, ...,λ λɶ ɶ . Thay các giá trị riêng này vào hệ thuần nhất tương ứng, ta được các hàm số riêng tương ứng độc lập tuyến tính của ( )K t,s . Ví dụ 2.1. Giải xấp xỉ phương trình tích phân Fredholm loại hai với nhân khơng suy biến ( ) ( ) ( )1 t s t 0 f t t 1 e f s ds e t= − + −∫ . Dùng cơng thức cầu phương Simpson ( ) ( ) ( ) ( ) 1 0 1y x dx y 0 4y 0.5 y 1 6 ≈  + +  ∫ , để tìm nghiệm gần đúng tại 1 2 3t 0 ; t 0.5 ; t 1= = = , ta được hệ phương trình tuyến tính (theo (2.6)) ( ) ( ) ( ) 11 1 12 2 13 3 1 21 1 22 2 23 3 2 31 1 32 2 33 3 3 1 K f K f K f g K f 1 K f K f g K f K f 1 K f g  − − − =   − + − − =  − − + − = ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ (2.7) Do t t s 1 2 3 1 2 1g(t) e t ;w ,w ,w ;K(t,s) t(1 e ) 6 3 6 = − = = = = − ta suy ra 19 ( ) ( ) ( ) ( ) 11 12 13 21 0,25 22 0,5 23 0,5 31 32 33 K K(0,0) 0;K K(0,0.5) 0;K K(0,1) 0 K K(0.5,0) 0; 1K K(0.5,0.5) 1 e ; 2 1K K(0.5,1) 1 e 2 K K(1,0) 0 ; K K(1,0.5) 1 e ; K K(1,1) 1 e  = = = = = =   = =   = = −   = = −  = = = = − = = − Từ đĩ, ta cĩ ( ) ( ) ( ) ( ) 11 12 13 0,25 0,5 21 22 23 0,5 31 32 33 K K K 0 1 1K 0 ; K 1 e ; K 1 e 3 12 2 1K 0 ; K 1 e ; K 1 e 3 6   = = =   = = − = −   = = − = − ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ và 0.5 1 2 3g 1 ; g e 0.5 ; g e 1= = − = − . Thay vào (2.7) ta được ( ) ( ) ( ) ( ) 1 0.25 0.5 0.5 2 3 0.5 2 3 f 1 1 1 e 2 f e 1 f e 0.5 3 12 2 1 e 1 f e 5 f e 1 3 6   =   + + − = −   − + + = − (2.8) Với các xấp xỉ 0.25 0.5e 1.2840254 ; e 1.6487213≈ ≈ ( ) ( ) ( ) ( ) 0.25 0.5 0.5 1 1 e 2 1.0947 ; e 1 0.0541 3 12 2 1 e 1 0.4325 ; e 5 1.2864 3 6 + ≈ − ≈ − ≈ + ≈ 20 (2.8) trở thành 1 2 3 2 3 f 1 1.0947f 0.0541f 1.1487 0.4325f 1.2864f 1.7183 =  + =  + = và ta được 1 2 3f 1 ; f 0.9999 ; f 0.9996= = = . Ta nhận được một biểu thức giải tích của nghiệm xấp xỉ theo cơng thức (2.3) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )1 1 1 2 2 2 3 3 3f t w K t,s f s w K t,s f s w K t,s f s g t= + + + trong đĩ 1 2 3 0.5 t 1 2 1 w , w , w 6 3 6 K(t,0) 0, K(t,0.5) t(1 e ) ,K(t,1) t(1 e ) f (0) 1, f (0.5) 0.9999, f (1) 0.9996  = = =  = = − = −  = = =   và do đĩ ( ) ( ) ( )0.5t t t2 1f t t 1 e 0.9999 t 1 e 0.9996 e t3 6= − + − + − hay ( ) 0.5t t tf t 0.1668t 0.6666te 0.1666te e= − − − + Nghiệm này cho các giá trị xấp xỉ tại một vài điểm như sau ( ) ( ) ( ) ( ) ( )f 0 1 ; f 0.5 1.0000 ; f 1 0.9996; f 0.4 1.0000 ; f 0.8 0.9999= = = = = . Chú ý rằng nghiệm chính xác của phương trình là ( )f t 1= . 21 1.2. Phương pháp thay nhân bằng nhân suy biến Khi nhân ( )K t,s ở phương trình (2.1) cĩ thể biểu diễn dưới dạng ( ) ( ) ( )n i i i 1 K t,s A t B s = =∑ (2.9) thì nhân ( )K t,s được gọi là nhân suy biến, chẳng hạn các nhân ( ) t st s, t s, sin t s , e ,...−+ ⋅ + là các nhân suy biến. Chú ý rằng nhân ( ) t sK t,s e= tuy là nhân khơng suy biến nhưng ta cĩ thể xấp xỉ ( )K t,s bằng nhân suy biến với độ chính xác tùy ý bằng khai triển Taylor ( )knt s k 0 t s e k! = =∑ . Giả sử ta phải giải phương trình (2.1) với nhân dạng (2.9). Khi đĩ, ta giả sử các hệ ( ) ( ) ( ){ }1 2 nA t ,A t ,...,A t và ( ) ( ) ( ){ }1 2 nB s ,B s ,...,B s là các hệ độc lập tuyến tính trên [ ]a ; b (vì trong trường hợp ngược lại, ta cĩ thể rút gọn bớt số hạng trong tổng ở vế phải của (2.9)). Thay (2.9) vào phương trình (2.1), ta được ( ) ( ) ( ) ( ) b n i i i 1a f t A (t)B s f s ds g t = = λ +∑∫ hay ( ) ( ) ( ) ( ) ( ) bn i i i 1 a f t A t B s f s ds g t = = λ +∑ ∫ (2.10) Đặt ( ) ( ) b i i a C B s f s ds= ∫ (2.11) thì nghiệm của phương trình đã cho cĩ thể được tìm dưới dạng 22 ( ) ( ) ( )n i i i 1 f t g t C A t = = + λ∑ (2.12) Thay (2.12) vào (2.10), ta cĩ ( ) ( ) ( ) ( ) ( ) ( ) ( ) bn n n i i i i j j i 1 i 1 j 1a g t C A t A t B s g s C A s ds g t = = =   + λ = λ + λ +    ∑ ∑ ∑∫ Rút gọn hai vế, ta suy ra ( ) ( ) ( ) ( ) ( ) bn n n i i i i j j i 1 i 1 j 1a C A t A t B s g s C A s ds = = =     = + λ       ∑ ∑ ∑∫ nghĩa là ( ) ( ) ( ) ( ) ( ) ( ) ( ) b bn n n n i i i i i j i j i 1 i 1 i 1 j 1a a C A t A t B s g s ds A t C B s A s ds = = = = = + λ∑ ∑ ∑ ∑∫ ∫ (2.13) Đặt ( ) ( ) b i i a g B s g s ds= ∫ , ( ) ( ) b ij i j a B s A s dsα = ∫ và vì ( ) ( ) ( ){ }1 2 nA t ,A t ,...,A t độc lập tuyến tính, (2.13) trở thành hệ phương trình tuyến tính n i i j ij j 1 C g C = = + λ α∑ hay n i ij j i j 1 C C g = − λ α =∑ (2.14) Nếu định thức của (2.14) là ( )D 0λ ≠ thì ta tìm được duy nhất các giá trị iC ; i 1,...,n= và ta cĩ nghiệm duy nhất của phương trình ở cơng thức (2.12). 23 Nếu ( )D 0λ = khi kλ = λ , thì kλ là giá trị riêng của ( )K t,s và ta tìm nghiệm của (2.14) ứng với ig 0= , ta sẽ được hệ hàm số riêng độc lập tuyến tính của ( )K t,s tương ứng với kλ . Như vậy, thuật tốn giải phương trình (2.1) với nhân suy biến gồm các bước sau : Bước 1. Tính các tích phân ( ) ( ) b i i a g B s g s ds= ∫ ; ( ) ( ) b ij i j a B s A s dsα = ∫ ; ( )i, j 1,...,n= Bước 2. Giải hệ phương trình tuyến tính n i ij j i j 1 C C g = − λ α =∑ ; ( )i 1,...,n= Nếu định thức của hệ khác khơng, thì tồn tại duy nhất iC ( )i 1,...,n= Bước 3. Nghiệm cần tìm cĩ dạng ( ) ( ) ( )n i i i 1 f t g t C A t = = + λ∑ . Ví dụ 2.2. Giải phương trình tích phân Fredholm loại hai với nhân suy biến, ( ) ( ) ( )1 t s t 0 f t e f s ds e t+= − + +∫ . Ta cĩ nhân ( ) t s t sK t,s e e .e+= − = − cho ( ) tA t e= − và ( ) sB s e= . Ta tìm nghiệm của phương trình dưới dạng, ( ) ( )t t tf t e t C A t e t Ce= + + ⋅ = + − . Ta cĩ ( ) ( )1 1 1s s s 2s 2 0 0 0 1g s e e ds se ds e ds e 1 2 = + = + = +∫ ∫ ∫ 24 ( )1 1s s 2s 2 0 0 1 e e ds e ds 1 e 2 α = − = − = −∫ ∫ và ta cĩ phương trình ( ) ( )2 21 1C e 1 C 1 e C 1 2 2 = + + − ⇔ = Vậy, nghiệm của phương trình là t tf (t) e t e t= + − = . Ví dụ 2.3. Ta khảo sát lại ví dụ 2.1 với phương trình ( ) ( ) ( )1 t s t 0 f t t 1 e f s ds e t= − + −∫ . Ta xấp xỉ ( )K t,s bằng cơng thức Taylor, 2 2 3 3 2 3 t s 2 3 4t s t s s st (1 e ) t 1 1 ts t ( s) t t 2 6 2 6        − ≈ − + + + = − + − + −              Thay vào phương trình ban đầu, ta được ( ) ( ) 1 2 3 2 3 4 t 0 s sf t t ( s) t t f s ds e t 2 6      = − + − + − + −          ∫ và ta nhận được phương trình Fredholm loại hai với nhân suy biến. Ta tìm nghiệm của phương trình dưới dạng ( ) t 2 3 41 2 3f t e t C t C t C t= − + + + , trong đĩ 1 2 3C ,C ,C là nghiệm của hệ (theo (2.14)) 1 1 11 2 12 3 13 1 2 1 21 2 22 3 23 2 3 1 31 2 32 3 33 3 C (C C C ) g C (C C C ) g C (C C C ) g − α + α + α =  − α + α + α =  − α + α + α = hay 25 11 1 12 2 13 3 1 21 1 22 2 23 3 2 31 1 32 2 33 3 3 (1 )C C C g C (1 )C C g C C (1 )C g − α − α − α =  −α + − α − α =  −α − α + − α = Ta được 1 s 1 0 1 2g (e s)sds 1 3 3   = − − = − − = −    ∫ 1 2 s 2 0 s e 9g (e s) ds 2 2 8 = − − = − +∫ 1 3 s 3 0 s 2e 29g (e s) ds 6 6 30 = − − = −∫ 1 1 1 3 4 5 11 12 13 0 0 0 1 1 1 s ds ; s ds ; s ds 4 5 6 α = − = − α = − = − α = − = −∫ ∫ ∫ 1 1 1 4 5 6 21 22 23 0 0 0 1 1 1 1 1 1 s ds ; s ds ; s ds 2 10 2 12 2 14 α = − = − α = − = − α = − = −∫ ∫ ∫ 1 1 1 5 6 7 31 32 33 0 0 0 1 1 1 1 1 1 s ds ; s ds ; s ds 6 36 6 42 6 48 α = − = − α = − = − α = − = −∫ ∫ ∫ Thay các ig và ijα vào hệ phương trình trên, ta được 1 2 3 1 2 3 1 2 3 5 1 1 2C C C 4 5 6 3 1 13 1 e 9C C C 10 12 14 2 8 1 1 49 2e 29C C C 36 42 48 6 30  + + = −   + + = − +   + + = −  Giải hệ này, ta được (lấy e 2.7183≈ ) 26 1 2 3 C 0.5010 C 0.1671 C 0.0418 = −  = −  = − và ta nhận được một biểu thức giải tích của nghiệm xấp xỉ ( ) t 2 3 4f t e t 0.5010t 0.1671t 0.0418t= − − − − . Các giá trị xấp xỉ của nghiệm tại một vài giá trị là ( ) ( ) ( ) ( ) ( )f 0 1 ; f 0.5 1.0000 ; f 1 1.0084; f 0.4 0.9999 ; f 0.8 1.0022= = = = = . Nhắc lại rằng, nghiệm chính xác của phương trình là ( )f t 1= . 2. CHƯƠNG TRÌNH TÍNH XẤP XỈ NGHIỆM PHƯƠNG TRÌNH TÍCH PHÂN FREDHOLM LOẠI HAI Trong phần này, chúng tơi đưa ra hai chương trình dùng để giải các phương trình Fredholm tuyến tính loại hai. Chương trình fred2 được cài đặt giải hệ (2.6) bằng phân rã LU. Chương trình fredin nhằm tính nội suy Nystrom của phương trình (2.3). 2.1. Cấu trúc chương trình Các chương trình được sử dụng gồm fred2 và fredin. Chương trình fred2 cĩ dữ liệu đầu vào là một con trỏ hàm g trả về ( )g t và một con trỏ hàm khác ( )ak trả về ( )K t,s , dữ liệu đầu ra là hàm nghiệm f tại các điểm quadrature, các điểm quadrature và trọng số. Các dữ liệu ra này được sử dụng cho chương trình fredin nhằm tính nội suy Nystrom của phương trình (2.3) từ đĩ trả về giá trị của f tại bất kỳ điểm nào trong đoạn [ ]a ; b . 2.1.1. Chương trình fred2 Chương trình này nhằm giải phương trình tích phân Fredholm loại hai cĩ dạng 27 ( ) ( ) ( ) ( ) b a f t K t,s f s ds g t= λ +∫ . Input : a, b : Giới hạn trên và giới hạn dưới của tích phân, n : Số điểm được sử dụng trong luật xấp xỉ tích phân, con trỏ hàm g : Trỏ đến hàm ( )g t , con trỏ hàm ak : Trỏ đến hàm ( )K t,s . Output : mảng t cĩ n phần tử : Các hệ số it của luật xấp xỉ tích phân Gaussian, mảng f cĩ n phần tử : Giá trị hàm nghiệm f tại các điểm it , mảng w cĩ n phần tử : Các trọng số Gaussian iw . Mã lệnh chương trình : void fred2(int n, float a, float b, float t[], float f[], float w[],float (*g)(float), float (*ak)(float, float)) { int i,j,*indx; (1) float d,**omk; indx=ivector(1,n); (2) omk=matrix(1,n,1,n); (3) gauleg(a,b,t,w,n); (4) for (i=1;i<=n;i++) { (5) for (j=1;j<=n;j++) (6) omk[i][j]=(float)(i == j)-(*ak)(t[i],t[j])*w[j]; f[i]=(*g)(t[i]); (7) } ludcmp(omk,n,indx,&d); (8) lubksb(omk,n,indx,f); (9) 28 free_matrix(omk,1,n,1,n); (10) free_ivector(indx,1,n); } Chú thích (1) : Khai báo các biến sử dụng trong chương trình. (2) : Gán indx là vector (mảng một chiều) cĩ n phần tử, indx được sử dụng cho bước (8). (3) : Gán omk là ma trận (mảng hai chiều) cĩ kích thước 1..n 1..n× . (4) : Gọi chương trình gauleg để lấy các hệ số it và trọng số iw của cơng thức xấp xỉ tích phân cho n điểm Gauss-Legendre. (5), (6) : Vịng lặp hai chiều để gán ma trận omk là ma trận ( )1 K− λ ɶ . (5), (7) : Với i chạy từ 1 đến n, gán [ ]f i bằng ( )ig t . Lúc này f sẽ là cột tự do g trong hệ ( )1 K f g− λ ⋅ =ɶ . (8) : Gọi chương trình ludcmp để thực hiện phân rã LU ma trận hệ số của hệ ( )1 K f g− λ ⋅ =ɶ . Chương trình này nhận vào ma trận vuơng omk cấp n, và trả về kết quả phân rã trở lại trong omk, mảng chứa chỉ số indx và d. Các kết quả trả về này sẽ được sử dụng cho chương trình lubksb để thực hiện giải hệ với một cột tự do cụ thể. (9) : Gọi chương trình lubksb để giải hệ với cột tự do chứa trong f (đã được gán ở bước (7)). Nghiệm trả về được lưu trở lại trong f. Lúc này, [ ]f i cho ta giá trị của hàm nghiệm f tại các điểm quadrature it . Các giá trị này sẽ là căn cứ cho chương trình fredin thực hiện nội suy Nystrom để tính giá trị hàm nghiệm f tại một điểm bất kỳ trong đoạn [ ]a ; b . (10) : Giải phĩng các biến trung gian, và kết thúc chương trình. Các kết quả output của chương trình lúc này được chứa trong các mảng t, f và w. 29 2.1.2. Chương trình fredin Chương trình này nhằm nội suy giá trị của hàm nghiệm f tại một điểm x bất kỳ trong đoạn [ ]a ; b bằng cơng thức nội suy Nystrom (2.3). Input : a, b : Giới hạn trên và giới hạn dưới của tích phân, n : Số điểm được sử dụng trong luật xấp xỉ tích phân, con trỏ hàm g : Trỏ đến hàm ( )g t , con trỏ hàm ak : Trỏ đến hàm ( )K t,s , mảng t cĩ n phần tử : Các hệ số it của luật xấp xỉ tích phân Gaussian, mảng f cĩ n phần tử : Giá trị hàm nghiệm f tại các điểm it , mảng w cĩ n phần tử : Các trọng số Gaussian iw . (t, f, và w là các Output của chương trình fred2) Output : Giá trị nội suy. Mã lệnh chương trình : float fredin(float x, int n, float a, float b, float t[], float f[], float w[], float (*g)(float), float (*ak)(float, float)) { int i; float sum=0.0; for (i=1;i<=n;i++) sum += (*ak)(x,t[i])*w[i]*f[i]; return (*g)(x)+sum; } 30 2.2. Kết quả thử nghiệm Ví dụ 2.4. Giải xấp xỉ phương trình tích phân Fredholm loại hai với nhân suy biến ( ) ( ) ( )1 t s t 0 f t e f s ds e t+= − + +∫ trong đĩ ( )f t là Nn hàm cần tìm, ( ) t sK t,s e += − là nhân, ( ) tg t e t= + . Phương trình này cĩ nghiệm chính xác là ( )f t t= . Mã lệnh chương trình DP g(const DP t) (1) { const DP EE = 2.71828; return (pow(EE,t)+t); } DP ak(const DP t, const DP s) (2) { const DP EE = 2.71828; return -pow(EE,t+s); } int main(void) { const int N=8; (3) DP a=0.0,b=1.0,x,ans; (4) Vec_DP t(N),f(N),w(N); (5) NR::fred2(a,b,t,f,w,g,ak); (6) cout << fixed << setprecision (6);(7) cout << "Vi du 1" << endl; for (;;) { (8) cout << "Nhap vao gia tri x trong khoang tu 0 den 1" << endl; (8.1) cin >> x; (8.2) if (x < 0.0) break; (8.3) 31 ans=NR::fredin(x,a,b,t,f,w,g,ak); (8.4) cout << "x da nhap, ket qua f(x) xap xi, ket qua f(x) thuc" << endl; cout << setw(10) << x << setw(11) << ans; (8.5) cout << setw(11) << x << endl << endl; (8.6) } return 0; (9) } Chú thích (1) : Định nghĩa hàm ( ) tg t e t= + . (2) : Định nghĩa hàm ( ) t sK t,s e += − . (3) : Khai báo hằng số N bằng 8, là số điểm của luật xấp xỉ tích phân. (4) : Khai báo các giới hạn trên, giới hạn dưới của tích phân, biến x, ans sẽ được dùng khi chạy thử nghiệm kết quả ở bước (8). (5) : Khai báo các vectơ t, f, w cĩ N phần tử. (6) : Gọi chương trình fred2, các kết quả trả về được lưu trong các vectơ t, f, w. (7), (8) : Yêu cầu người sử dụng nhập vào một giá trị x để thử nghiệm (8.1, 8.2) (nếu x bé hơn 0 thì kết thúc chương trình (8.3)). Sau đĩ chương trình sẽ gọi chương trình fredin tính nội suy Nystrom để cĩ xấp xỉ của ( )f x (được trả về trong biến ans) (8.4), và xuất các kết quả gồm : giá trị x được nhập vào, giá trị xấp xỉ ( )f x (8.5), giá trị ( )f x thực (trong ví dụ này, hàm nghiệm chính xác là ( )f x x= ) (8.6). (9) : Kết thúc chương trình. Với một số giá trị x = 0.2 ; 0.5 ; 0.8, chương trình nêu trên cho ta kết quả như sau : 32 Ví dụ 2.5. Giải xấp xỉ phương trình tích phân Fredholm loại hai với nhân khơng suy biến ( ) ( ) ( )1 t s t 0 f t t 1 e f s ds e t= − + −∫ trong đĩ, ( )f t là Nn hàm cần tìm, ( ) ( )t sK t,s t 1 e= − là nhân, ( ) tg t e t= − . Phương trình này cĩ nghiệm chính xác là ( )f t 1= . Mã lệnh chương trình DP g(const DP t) { const DP EE = 2.71828; return (pow(EE,t)-t); } DP ak(const DP t, const DP s) { const DP EE = 2.71828; return t*(1-pow(EE,t*s)); } int main(void) { const int N=8; DP a=0.0,b=1.0,x,ans; Vec_DP t(N),f(N),w(N); NR::fred2(a,b,t,f,w,g,ak); 33 cout << fixed << setprecision(6); cout << "Vi du 2" << endl; for (;;) { cout << "Nhap vao gia tri x trong khoang tu 0 den 1" << endl; cin >> x; if (x < 0.0) break; ans=NR::fredin(x,a,b,t,f,w,g,ak); cout << "x da nhap, ket qua f(x) xap xi, ket qua f(x) thuc" << endl; cout << setw(10) << x << setw(11) << ans; cout << setw(11) << 1 << endl << endl; } return 0; } Với một số giá trị x = 0.2 ; 0.4 ; 0.7, chương trình cho kết quả sau Ví dụ 2.6. Giải xấp xỉ phương trình tích phân Fredholm loại hai với nhân dạng tích chập ( ) ( ) 1 t t s 0 f t e t e f s ds−= + − ∫ trong đĩ, ( )f t là Nn hàm cần tìm, ( ) ( )t sK t,s e t s−= − = ϕ − là nhân, ( ) tg t e t= + . Phương trình này cĩ nghiệm chính xác là ( ) t 1f t t e −= + . 34 Mã lệnh chương trình DP g(const DP t) { const DP EE = 2.71828; return (pow(EE,t)+t); } DP ak(const DP t, const DP s) { const DP EE = 2.71828; return -pow(EE,t-s); } int main(void) { const DP EE = 2.71828; const int N=8; DP a=0.0,b=1.0,x,ans; Vec_DP t(N),f(N),w(N); NR::fred2(a,b,t,f,w,g,ak); cout << fixed << setprecision(6); cout << "Vi du 3" << endl; for (;;) { cout << "Nhap vao gia tri x trong khoang tu 0 den 1" << endl; cin >> x; if (x < 0.0) break; ans=NR::fredin(x,a,b,t,f,w,g,ak); cout << "x da nhap, ket qua f(x) xap xi, ket qua f(x) thuc" << endl; cout << setw(10) << x << setw(11) << ans; cout << setw(11) << (x+(1.0/EE)*pow(EE,x)) << endl << endl; } return 0; } Với một số giá trị x = 0.2 ; 0.4 ; 0.7 ; 0.9, chương trình cho kết quả 35 Ví dụ 2.7. Giải xấp xỉ phương trình tích phân Fredholm loại hai với nhân dạng tích chập và cĩ nhiễu tham số ( ) ( ) ( ) 1 2 0 1 f t t s f s ds t t n = − + +∫ hay ( ) ( ) ( ) ( )1 2 0 f t n t s f s ds n t t= − + +∫ , trong đĩ ( )f t là Nn hàm cần tìm, ( ) ( ) ( )K t,s n t s t s= − = ϕ − là nhân, ( ) ( )2g t n t t= + . Phương trình này cĩ nghiệm chính xác là ( ) ( ) 3 2 3 2 2 2 2 n 10n 12n n 42nf t nt t n 12 6 n 12  − + +  − = + + + +  . Với n 10= ta cĩ ( ) 2 15 200f t 10t t 14 42 = + − . Với n 100= ta cĩ ( ) ( )2f t 100t 89,77227327 t 9,65508057= − + . Với n 1000= ta cĩ ( ) ( )2 957,9885041f t 1000t 989,9761203 t 6 = − + . 36 Mã lệnh chương trình DP n = 100; DP g(const DP t) { return n*(t*t+t); } DP ak(const DP t, const DP s) { return n*(t-s); } int main(void) { const int N=8; DP a=0.0,b=1.0,x,ans; Vec_DP t(N),f(N),w(N); NR::fred2(a,b,t,f,w,g,ak); cout << fixed << setprecision(6); cout << "Vi du 4 voi N = " << (int)n << endl; for (;;) { cout << "Nhap vao gia tri x trong khoang tu 0 den 1" << endl; cin >> x; if (x < 0.0) break; ans=NR::fredin(x,a,b,t,f,w,g,ak); cout << "x da nhap, ket qua f(x) xap xi, ket qua f(x) thuc" << endl; cout << setw(10) << x << setw(11) << ans; DP tmp = n*x*x + ((-n*n*n+10*n*n+12*n)/(n*n+12))*x + (n*n*n- 42*n*n)/(6*(n*n+12)); cout << setw(11) << tmp << endl << endl; } return 0; } Với n 10= và các giá trị x = 0.2 ; 0.4 ; 0.7 ; 0.9, chương trình cho 37 Với n 100= và các giá trị x = 0.2 ; 0.4 ; 0.7 ; 0.9, ta được Với n 1000= và các giá trị x = 0.2 ; 0.4 ; 0.7 ; 0.9, ta cĩ 38 Chương 3 THUẬT TỐN TOEPLITZ Trong chương này, chúng ta trình bày giải thuật giải số hệ phương trình tuyến tính AX B= , với A là ma trận Toeplitz đối xứng nhằm ứng dụng vào việc giải số các phương trình Fredholm loại một và Fredholm loại hai cĩ nhân là các hàm dạng tích chập. 1. THUẬT TỐN TOEPLITZ Trước hết, ta xét hệ phương trình 0 0 1 1 n n 0 1 0 0 1 n 1 n 1 n 0 n 1 1 0 n n t f t f t f g t f t f t f g t f t f t f g − − + + + =  + + + =    + + + = … … … … (3.1) hay dưới dạng ma trận nT f g= với ma trận 0 1 n 1 0 n 1 n n n 1 0 t t t t t t T t t t − −       =       … … … … là ma trận Toeplitz đối xứng cấp n 1+ và ma trận cột [ ]T0 1 ng g g ... g= là những ma trận đã biết, cịn ma trận cột [ ]T0 1 nf f f ... f= là ma trận Nn hàm cần tìm. Ta xây dựng [ ]Tn n0 n1 nnf f f f= … là nghiệm của hệ phương trình (3.1) bằng quy nạp theo n như sau 39 1. Với k 0= , đặt 00 0 0 0 1 0 00 0 00 1 0 a 1, t , t gf , f t t = α = β =   = γ =  2. Giả sử ta cĩ [ ]Tk k0 k1 kkf f f f= … là nghiệm của hệ (3.1) với n k= , nghĩa là ta cĩ k0 k1 kk k ka , a , , a ; ,α β… và k0 k1 kk kf , f , , f ; γ… sao cho các hệ (3.2) và (3.3) sau đây thỏa k0 0 k1 1 kk k k k0 1 k1 0 kk k 1 k0 k k1 k 1 kk 0 k0 k 1 k1 k kk 1 k a t a t a t a t a t a t 0 a t a t a t 0 a t a t a t − − + + + + = α  + + + =   + + + =  + + + = β … … … … … (3.2) k0 0 k1 1 kk k 0 k0 1 k1 0 kk k 1 1 k0 k k1 k 1 kk 0 k k0 k 1 k1 k kk 1 k f t f t f t g f t f t f t g f t f t f t g f t f t f t − − + + + + =  + + + =   + + + =  + + + = γ … … … … … (3.3) Ta cần chứng minh T k 1 k 1,0 k 1,1 k 1,k 1f f f f+ + + + + =  … là nghiệm của (3.1) với n k 1= + . a) Trước hết, ta tính k 1,0 k 1,1 k 1,k k 1,k 1 k 1 k 1a , a , , a , a ; ,+ + + + + + +α β… bằng cách đặt 40 k 1,0 k0 k 1,1 k1 k kk k 1,k kk k k1 k 1,k 1 k k0 a a a a p a a a p a a p a + + + + + =  = +    = +   = … (3.4) và k 1 k k k k 1 k 1,0 k 2 k 1,1 k 1 k 1,k 1 1 p a t a t a t + + + + + + + + α = α + β β = + + + … (3.5) với k k k p β= − α (3.6) thì ta cần chứng tỏ rằng k 1,0 0 k 1,1 1 k 1,k 1 k 1 k 1 k 1,0 1 k 1,1 0 k 1,k 1 k k 1,0 k 1 k 1,1 k k 1,k 1 0 k 1,0 k 2 k 1,1 k 1 k 2,k 2 1 k 1 a t a t a t a t a t a t 0 a t a t a t 0 a t a t a t + + + + + + + + + + + + + + + + + + + + + + + + + = α  + + + =    + + + =   + + + = β … … … … … (3.7) Thật vậy, để kiểm chứng (3.7), ta cĩ k 1,0 0 k 1,1 1 k 1,k k k 1,k 1 k 1a t a t a t a t+ + + + + ++ + + +… k0 0 k1 k kk 1 kk k k1 k k k0 k 1a t (a p a )t (a p a )t p a t += + + + + + +… k0 0 k1 1 kk k k k0 k 1 k1 k kk 1(a t a t a t ) p (a t a t a t )+= + + + + + + +… … k k k k 1p += α + β = α k 1,0 1 k 1,1 0 k 1,k k 1 k 1,k 1 ka t a t a t a t+ + + − + ++ + + +… k0 1 k1 k kk 0 kk k k1 k 1 k k0 ka t (a p a )t (a p a )t p a t−= + + + + + +… 41 k0 1 k1 0 kk k 1 k k0 k k1 k 1 kk 0(a t a t a t ) p (a t a t a t )− −= + + + + + + +… … k0 p 0 0= + = k 1,0 k 1 k 1,1 k k 1,k 1 k 1,k 1 0a t a t a t a t+ + + + + ++ + + +… k0 k 1 k1 k kk k kk k k1 1 k k0 0a t (a p a )t (a p a )t p a t+= + + + + + +… k0 k 1 k1 k kk 1 k k0 0 k1 1 kk k(a t a t a t ) p (a t a t a t )+= + + + + + + +… … k k k k 0β= β − α = α k 1,0 k 2 k 1,1 k 1 k 2,k 2 1 k 1a t a t a t+ + + + + + ++ + + = β… theo (3.5) b) Tiếp theo, ta tính k 1,0 k 1,1 k 1,k k 1,k 1 k 1f , f , , f , f ;+ + + + + +γ… bằng cách đặt k 1,0 k0 k k 1,k 1 k 1,1 k1 k k 1,k k 1,k kk k k 1,1 k 1,k 1 k k 1,0 f f q a f f q a f f q a f q a + + + + + + + + + + = +  = +    = +   = … (3.8) và k 1 k 1,0 k 2 k 1,1 k 1 k 1,k 1 1f t f t f t+ + + + + + +γ = + + +… với k 1 k k k 1 gq + + − γ = α (3.9) Khi đĩ, ta cần phải chứng tỏ rằng 42 k 1,0 0 k 1,1 1 k 1,k 1 k 1 0 k 1,0 1 k 1,1 0 k 1,k 1 k 1 k 1,0 k 1 k 1,1 k k 1,k 1 0 k 1 k 1,0 k 2 k 1,1 k 1 k 1,k 1 1 k 1 f t f t f t g f t f t f t g f t f t f t g f t f t f t + + + + + + + + + + + + + + + + + + + + + + + + + =  + + + =    + + + =   + + + = γ … … … … … (3.10) Thật vậy, ta cĩ k 1,0 0 k 1,1 1 k 1,k k k 1,k 1 k 1f t f t f t f t+ + + + + ++ + + +… k0 k k 1,k 1 0 k1 k k 1,k 1 kk k k 1,1 k k k 1,0 k 1(f q a )t (f q a )t (f q a )t q a t+ + + + + += + + + + + + +… k0 0 k1 1 kk k k k 1,0 k 1 k 1,1 k k 1,k 1 0(f t f t f t ) q (a t a t a t )+ + + + += + + + + + + +… … 0 k 0g q .0 g= + = k 1,0 1 k 1,1 0 k 1,k k 1 k 1,k 1 kf t f t f t f t+ + + − + ++ + + +… k0 k k 1,k 1 1 k1 k k 1,k 0 kk k k 1,1 k 1 k k 1,0 k(f q a )t (f q a )t (f q a )t q a t+ + + + − += + + + + + + +… k0 1 k1 0 kk k 1 k k 1,0 k k 1,1 k 1 k 1,k 1 1(f t f t f t ) q (a t a t a t )− + + − + += + + + + + + +… … 1 k 1g q .0 g= + = k 1,0 k 1 k 1,1 k k 1,k 1 k 1,k 1 0f t f t f t f t+ + + + + ++ + + +… k0 k 1 k1 k kk 1 k k 1,0 0 k 1,1 1 k 1,k 1 k 1(f t f t f t ) q (a t a t a t )+ + + + + += + + + + + + +… … k 1 k k k k 1 k k 1 k 1 k 1 gq g++ + + +   − γ = γ + α = γ + α =  α  Từ đĩ, ta xây dựng được thuật tốn giải hệ phương trình (3.1) gồm các bước như sau : Bước 1. Với k = 0, đặt 43 00 0 0 0 1 0 00 0 00 1 0 a 1, t , t gf , f t t = α = β =   = γ =  Bước 2. Giả sử, ta đã xây dựng được k0 k1 kk k k k0 k1 kk k a , a , , a ; , f , f , , f ; α β  γ … … thì ta đặt k 1,0 k0 k 1,1 k1 k kk k 1,k kk k k1 k 1,k 1 k k0 a a a a p a a a p a a p a + + + + + =  = +    = +   = … k 1 k k kp+α = α + β với k k k k 1 k k k 1 p gq + + β = − α  − γ =  α Bước 3. Nghiệm của hệ phương trình là T k 1 k 1,0 k 1,1 k 1,k 1f f f f+ + + + + =  … với k 1,0 k0 k k 1,k 1 k 1,1 k1 k k 1,k k 1,k kk k k 1,1 k 1,k 1 k k 1,0 f f q a f f q a f f q a f q a + + + + + + + + + + = +  = +    = +   = … 44 2. CHƯƠNG TRÌNH GIẢI HỆ PHƯƠNG TRÌNH TOEPLITZ Chương trình mytoeplitz Input : n : Cấp của ma trận Toeplitz. t : Mảng một chiều gồm n phần tử it ( )i 0,...,n 1= − đại diện cho n đường chéo của ma trận Toeplitz đối xứng. g : Mảng một chiều gồm n phần tử ứng với vectơ cột hệ số tự do của hệ cần giải. Output : Mảng f cĩ n phần tử ứng với vectơ nghiệm tìm được. Mã lệnh chương trình : void NR::mytoeplitz(float[] t, float[] g, float[] f, int n) { int i; float alpha,beta,gamma,p,q; (1) float *a, *anext; (2) a = fvector(1,n); anext = fvector(1,n); if (n<2) (3) { nrerror("Size of n must be larger than 1"); return; } a[0] = 1; (4) alpha = t[0]; beta = t[1]; f[0] = g[0]/t[0]; gamma = f[0]*t[1]; for(int k=1;k<n;k++) (5) 45 { p = -beta/alpha; (5.1) anext[0] = a[0]; (5.2) for(i=1;i<k;i++) { anext[i] = a[i] + p*a[k-i]; } anext[k] = a[0]*p; alpha = alpha + p*beta; (5.3) beta = 0; (5.4) for(i=0;i<=k;i++) { beta += anext[i]*t[k+1-i]; } q = (g[k] - gamma)/alpha; (5.5) for(i=0;i<k;i++) (5.6) { f[i] = f[i] + q*anext[k-i]; } f[k] = q*anext[0]; gamma = 0; (5.7) for(i=0;i<=k;i++) { gamma += f[i]*t[k+1-i]; } for(i=0;i<=k;i++) { a[i] = anext[i]; } } } 46 Chú thích (1) : Khai báo các biến sử dụng trong chương trình. alpha ứng với kα , beta ứng với kβ , gamma ứng với kγ , p ứng với kp , q ứng với kq . (2) : Khai báo mảng a và anext cĩ n phần tử, [ ]ia (phần tử thứ i của mảng a, i 0..n 1= − ) ứng với kia , [i]anext ứng với k 1, ia + . (3) : Kiểm tra nếu n 2< thì thốt khỏi chương trình. (4) : Gán các giá trị khởi tạo : 00a 1= , 0tα = , 1tβ = , 000 0 gf t = , 0 00 1f tγ = . (5) : Cho vịng lặp với k chạy từ 1 đến n 1− , lần lượt tính các kia , kα , kβ , kif , kγ . (5.1) : Tính kk k p β= − α . (5.2) : Tính k 1,0 k0a a+ = , k 1,1 k1 k kka a p a+ = + , ..., k 1, k kk k k1a a p a+ = + , k 1,k 1 k0 ka a p+ + = . (5.3) : Tính k 1 k k kp+α = α + β . (5.4) : Tính k 1+β theo cơng thức k 1 k 1,0 k 2 k 1,1 k 1 k 1,k 1 1a t a t ... a t+ + + + + + +β = + + + (giá trị của k 1+β được gán trở lại vào biến beta vì phần sau ta khơng cịn dùng đến kβ nữa). (5.5) : Tính k 1 kk k 1 gq + + − γ = α . (5.6) : Tính các k 1,if + theo cơng thức : k 1,0 k0 k k 1,k 1f f q a+ + += + , k 1,1 k1 k k 1, kf f q a+ += + , ..., k 1,k kk k k 1,1f f q a+ += + , k 1, k 1 k k 1,0f q a+ + += . (5.7) : Tính k 1 k 1,0 k 2 k 1,1 k 1 k 1, k 1 1f t f t ... f t+ + + + + + +γ = + + + . 47 3. KẾT QUẢ THỬ NGHIỆM 3.1. Kết quả thử nghiệm đối với hệ Toeplitz Ví dụ 3.1. Xét lại hệ phương trình (3.1) nT f g= trong đĩ ma trận 0 1 n 1 0 n 1 n n n 1 0 t t t t t t T t t t − −       =       … … … … là ma trận Toeplitz đối xứng cấp n 1+ và ma trận cột [ ]T0 1 ng g g ... g= . a) Với n = 9 và ( )0 1 9 0.000001,0.111111,0.222222,0.333333,0.444444,t , t ,..., t 0.555556,0.666667,0.777778,0.888889,1.0   =     [ ]T0 1 ng g g ... g= T0.132120, 0.102953, 0.084845, 0.076632, 0.077274, 0.085841, 0.101500, 0.123500, 0.151179, 0.183940   =     Nếu giải bằng thuật giải Toeplitz thì nghiệm của hệ là T0.026777, 0.049764, 0.044531, 0.039848, 0.035658, f 0.031908, 0.028552, 0.025550, 0.022863, 0.010608   =     Nếu ta dùng phần mềm Maple 12 thì kết quả là (làm trịn 6 chữ số thập phân) T0.026778, 0.049766, 0.044528, 0.039848, 0.035662, f 0.031914, 0.028535, 0.025556, 0.022869, 0.010605   =     b) Với n = 9 và 48 ( )0 1 9 1.00, 1111.11, 2222.22, 3333.33, 4444.44,t , t ,..., t 5555.56, 6666.67, 7777.78, 8888.89, 10000.00 − − − −  =   − − − − −  [ ]T0 1 ng g g ... g= . Ta xét hai trường hợp sau : - Trường hợp 1 : với T1321.2, 1029.53, 848.45, 766.32, 772.74, g 858.41, 1014.98, 1235.0, 1511.79, 1839.4 − − − − −  =   − − − − −  nếu giải bằng thuật giải Toeplitz, nghiệm thu được là T0.026789, 0.049752, 0.044531, 0.039848, 0.035658, f 0.031908, 0.028553, 0.025550, 0.022859, 0.010615   =     và nếu ta dùng phần mềm Maple 12 thì kết quả là (làm trịn 6 chữ số thập phân) T0.026791, 0.049753, 0.044528, 0.039848, 0.035662, f 0.031905, 0.028553, 0.025550, 0.022864, 0.010612   =     - Trường hợp 2 : với [ ]Tg 1000, 2000,3000, 4000,5000, 6000,7000, 8000,9000, 10000= − − − − − nếu giải bằng thuật giải Toeplitz, nghiệm thu được là T1.79726, 3.593531, 5.390298, 7.187063, 8.983829, f 10.780595, 12.577361, 14.374125, 16.175334, 8.089396 − −  =   − − −  và nếu ta dùng phần mềm Maple 12 thì kết quả là (làm trịn 6 chữ số thập phân) T1.797212, 3.593430, 5.390182, 7.186933, 8.983721, f 10.780545, 12.577333, 14.374082, 16.175277, 8.089364 − −  =   − − −  49 3.2. Giải một số phương trình Fredholm bằng cách đưa về hệ Toeplitz Ví dụ 3.2. Cho phương trình Fredholm loại một ( ) 1 t 0 e 1 e 2 t s f s ds e t 2e 2e − + +  − = + −    ∫ (3.11) trong đĩ ( )f t là Nn hàm cần tìm, hàm nhân ( ) ( )K t,s t s t s= ϕ − = − và hàm t e 1 e 2g(t) e t 2e 2e − + +  = + −    . Nghiệm chính xác của phương trình là t1f (t) e 2 − = . a) Giải trực tiếp sử dụng thuật giải Toeplitz : Ta chia đoạn [ ]0; 1 thành n 1− đoạn bằng nhau [ ]0 1s ,s , [ ]1 2s ,s , ... , [ ]n 2 n 1s ,s− − , với 0s 0= , n 1s 1− = , nghĩa là mỗi đoạn cĩ độ dài 1h n 1 = − . Xấp xỉ tích phân ở vế trái của (3.11) bằng cơng thức cầu phương, ( ) ( ) 1 n 1 i i i 00 f s ds w f s − = =∑∫ , và để đơn giản trong việc tính tốn ta chọn iw 1= , với mọi i 0,n 1= − . Khi đĩ (3.11) trở thành ( ) ( ) ( )n 1 i i i 0 g t K t,s f s − = =∑ Tại điểm jt t= ; j 0,...,n 1= − , ta cĩ n 1 j j i i i 0 g(t ) K(t ,s )f (s ) − = =∑ (3.12) hay 50 ( ) ( ) ( ) ( ) ( )n 1 n 1j j i i j i i i 0 i 0 g s K s ,s f s s s f s − − = = = = ϕ −∑ ∑ Đặt ( )ji j i j it s s t −= ϕ − = , j jg(s ) g= , i if (s ) f= (3.12) cĩ thể viết lại dạng n 1 j j i i i 0 g t f − − = =∑ hay dưới dạng phương trình ma trận nT .f g= trong đĩ 0 1 2 (n 1) 1 0 1 (n 2) n 2 1 0 (n 3) n 1 n 2 n 3 0 t t t t t t t t T t t t t t t t t − − − − − − − − − − − −         =          … … … ⋮ ⋮ ⋮ ⋱ ⋮ … là ma trận Toeplitz cấp n và 0 1 n 1 g g g g −       =       ⋮ là những ma trận đã biết và 0 1 n 1 f f f f −       =       ⋮ là ma trận Nn phải tìm. 51 Do ( )t s t sϕ − = − là hàm đối xứng, nên nT là ma trận Toeplitz đối xứng. Khi đĩ, nT cĩ dạng 0 1 n 1 1 0 n 2 n n 1 n 2 0 t t ... t t t ... t T ... ... ... ... t t ... t − − − −       =       . Giải ví dụ (3.2) với n = 10. Khi đĩ, h 0.111111= is 0 h i h i 0.111111 i= + ⋅ = ⋅ = ⋅ . ( ) ( ) ( ) T0 1 9 T g g s ,g s ,...,g s 0.13212,0.102953,0.084845,0.076632,0.077274, 0.085841,0.1015,0.1235,0.151179,0.18394  =     =     ( )0 1 9 0.000001,0.111111,0.222222,0.333333,0.444444,t , t ,..., t 0.555556,0.666667,0.777778,0.888889,1.0   =     (ta đã cộng thêm 0.000001 vào 0t để phù hợp với thuật giải Toeplitz). Kết quả nghiệm f là T0.026777,0.049764,0.044531,0.039848,0.035658, f 0.031908,0.028552,0.02555,0.022863,0.010608   =     , với nghiệm chính xác của hệ là t1f (t) e 2 − = . Ta cĩ bảng so sánh giữa nghiệm f thực và nghiệm f giải được theo cách trên : 52 (Với n càng lớn, thì kết quả càng xấu đi) b) Giải bằng cách đưa về dạng Fredholm loại hai : Ta viết lại (3.11) thành ( ) ( ) 1 0 0 t s f s ds g t= − −∫ Cộng thêm ( )1 f t M vào vế trái (với M khá lớn) và đặt ( ) t e 1 e 2g t e t 2e 2e −  + +  = − + −      Ta cĩ ( ) ( ) ( ) 1 0 1 f t t s f s ds g t M = − +∫ (3.13) Ta coi như (3.13) xấp xỉ tương đương với (3.11) khi M khá lớn, nghĩa là nghiệm của (3.13) cũng là nghiệm của (3.11) và ngược lại. Giải (3.13) với chương trình fred2 đã trình bày ở trên, với M = 100, số khoảng chia là 10, sử dụng chương trình gauleg để xấp xỉ tích phân, khơng sử dụng thuật giải Toeplitz (đương nhiên, vì các khoảng chia khơng đều nhau), ta cĩ kết quả 53 (cột Abscissa là cột các it , cột Calc soln là cột nghiệm f đã tính được, cột True soln là cột f tính từ nghiệm thực t1f (t) e 2 − = ) Thực hiện lại với M 1000= , ta cĩ Với M 10000= , ta cĩ Nếu khơng dùng chương trình gauleg để chọn khoảng it và iw , mà chọn it cách đều nhau, cho tất cả iw 1= , thì kết quả sẽ rất xấu : Với M 10000= 54 Với M 10000= và số khoảng chia là 100 Với số khoảng chia lớn hơn nữa thì chương trình chạy rất chậm, do chương trình fred2 này vẫn đang sử dụng kiểu giải ma trận bằng phân rã LU. Trong trường hợp ví dụ này, ta vẫn cĩ thể sử dụng thuật giải Toeplitz, vì trong phương trình Fredholm loại hai, ma trận hệ số là ( )i j j1 K w− λ (theo phương trình (2.6)), nên khi các khoảng chia đều nhau, jw 1= với mọi j, và ( ) ( )K t,s t s t s= ϕ − = − là hàm đối xứng, thì ma trận hệ số này cũng cĩ dạng Toeplitz đối xứng, thậm chí cịn tốt hơn trường hợp Fredholm loại một, vì 0t đã khác 0 ( 0t 1= ). Với M 10000= và số khoảng chia là 100 55 Nhận thấy rằng kết quả khơng hề khác so với khi sử dụng kiểu giải ma trận bình thường (điều này cho thấy thuật giải Toeplitz của ta là chính xác). Tuy nhiên, ta đã cĩ thể chạy với số khoảng chia lớn hơn mà chương trình khơng bị chậm. Với M 10000= và số khoảng chia là 1000 Tuy nhiên, vấn đề cần hướng đến là xấp xỉ được nghiệm thực thì ta vẫn chưa đạt được. 56 KẾT LUẬN Luận văn đã trình bày được giải thuật cơ bản cho việc xấp xỉ nghiệm của phương trình tích phân với dữ liệu rời rạc bằng việc khảo sát một hệ phương trình tuyến tính. Đối với phương trình tích phân Fredholm loại hai, ta được bài tốn chỉnh nhưng nếu cần xấp xỉ nghiệm tại một số lớn các điểm nút, hệ phương trình tuyến tính nhận được là tương đối lớn và như vậy việc sử dụng các giải thuật thơng thường như giải thuật LU sẽ cho thời gian thực khi chạy chương trình là chậm, thậm chí khơng khả thi. Luận văn khảo sát một loại phương trình tích phân quan trọng trong khá nhiều ứng dụng với nhân cĩ dạng ( ) ( )K t,s t s= ϕ − , với hàm số thực ϕ cho trước. Phương trình tích phân với nhân loại này cho hệ phương trình tuyến tính với ma trận các hệ số thuộc loại ma trận Toeplitz mà luận văn đưa ra giải thuật giải với độ phức tạp ( )O n thay vì ( )2O n cho trường hợp tổng quát. Đối với phương trình tích phân loại một, luận văn cũng đã đưa ra các tính số nghiệm xấp xỉ tương ứng với nghiệm của một số các phương trình tích phân loại hai. Kết quả nhận được của luận văn chưa thật sự thuyết phục do các phép xấp xỉ tích phân với các trọng số iw 1= và các điểm nút đều nhau khơng cho được sai số nhỏ cần thiết. Hướng phát triển bổ sung là các kết quả lý thuyết cho việc đánh giá các sai số trong việc xấp xỉ phương trình tích phân bằng hệ phương trình tuyến tính cũng như khảo sát trường hợp xấp xỉ tích phân bằng các trọng số cũng như các điểm nút thay đổi. 57 TÀI LIỆU THAM KHẢO [1] Đặng Đình Áng, Trần Lưu Cường, Huỳnh Bá Lân, Nguyễn Văn Nhân (2001), Biến Đổi Tích Phân, NXB Giáo Dục. [2] Phạm Kỳ Anh (2000), Giải Tích Số, NXB ĐHQG Hà Nội. [3] Nguyễn Minh Chương, Nguyễn Văn Khải, Khuất Văn Ninh, Nguyễn Văn Tuấn, Nguyễn Tường (2002), Giải Tích Số, NXB Giáo Dục. [4] Nguyễn Cơng Tâm, Đinh Ngọc Thanh (2005), Phương Trình Tích Phân, NXB ĐHQG TP. Hồ Chí Minh. [5] Lê Đình Thịnh (1995), Phương Pháp Tính, NXB Khoa học và Kỹ thuật Hà Nội. [6] Hồ Thuần, Nguyễn Cơng Thúy, Nguyễn Quý Hỷ, Phan Văn Hạp (1970), Cơ Sở Phương Pháp Tính, NXB Đại học và Trung học Chuyên nghiệp. [7] Trần Văn Trân (2007), Phương pháp số thực hành, NXB ĐHQG Hà Nội. [8] Dương Thủy Vỹ (2001), Phương Pháp Tính, NXB Khoa học và Kỹ thuật Hà Nội. [9] Andreas Kirsch (1996), An Introduction to the Mathematical Theory of Inverse Problems. [10] Andrei D. Polyanin and Alexander V. Manzhirov (1998), Handbook of Integral Equations, CRC Press Boca Raton London New York Washington. [11] E. C. Titchmarsh (1975), Introduction to the theory of Fourier integrals, 7th edition, Oxford University Press, London. [12] N.S. Bakhvalov (1977), Numerical Methods, Mir Publishers, Moscow. [13] Robert M. Gray, Toeplitz and Circulant Matrices: A review, Stanford University, Stanford 94305, USA. 58 [14] Seymour Lipschutz, Marc Lars Lipson (2001), Theory and Problems of Linear Algebra, Third edition, McGraw-Hill. [15] U. Grenander, G. Szegư (1958), Toeplitz Forms and Their Applications, Univ. of California Press, Berkeley and Los Angeles, 3, pp. 64-65. [16] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery (1992), Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press.