Luận văn Phương trình tích chập và hệ Toeplitz

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 hoàn toà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 hoà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 toá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 toá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 hoá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 toán sau : Bài toá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 toá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 toá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 toá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 TOÁ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 TOÁ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 toá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ì thoá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 toá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 toá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.