Học Spline với biên ràng buộc và biên tự nhiên

HỌC SPLINE VỚI BIÊN RÀNG BUỘC VÀ BIÊN TỰ NHIÊN NGUYỄN PHÚ VINH (*) TÓM TẮT Trong bài báo này, đầu tiên chúng tôi khảo sát họ đường cong spline, có ứng dụng nhiều trong toán học, cơ học, v.v.. Phương trình các đường cong spline là bậc 3, thường dùng để xấp xỉ hàm, như phương pháp phần tử hữu hạn thường sử dụng. Trong phương pháp này, chúng tôi thiết lập thuật toán, trưng ra cách xử lí điều kiện biên tự nhiên và điều kiện biên ràng buộc. Sau đó dùng matlab, thiết lập chương trình, với dữ liệu đưa vào các mốc nội suy rất đơn giản để tự động hoá đưa ra các biểu thức hình thức (symbolic) của các phương trình, và biểu diễn bằng đồ thị, phản ánh đúng với điều kiện biên cho trước. ABSTRACT In this paper, we first survey their spline curves which has many applications in mathematics, mechanical, etc. The equation of the spline curves is order 3 that are often used to appoximate the function such as finite element method. In this method, we set the algorithm, and present a way to handle boundary conditions binding. Then, the use of Matlab, and the setting of the program, with data interpolated into markets, are simple to make automated form of expressions and programs represented by graphs which reflect conditions the previous record. 1. MỞ ĐẦU Trong các phương pháp toán xấp xỉ hàm (ví dụ phương pháp phần tử hữu hạn, sai phân hữu hạn, v.v.), ta cần xấp xỉ hàm phải đi qua tất cả các mốc nội suy, và hàm xấp xỉ phải khả vi liên tục đến cấp 2. Đã có các phương pháp xấp xỉ Lagrange, Newton, nhưng đối với xấp xỉ nhiều (n) điểm mốc nội suy thì hai phương pháp này không còn thích hợp, vì chúng phải dùng đa thức cấp n-1 (bậc rất cao) để xấp xỉ, vì thế sai số sẽ lớn. Ngày nay, để khắc phục điều trên, với công nghệ tính toán hiện đại, người ta đã thay thế chúng bằng các hàm spline phức tạp hơn trong tính toán, nhưng xấp xỉ tốt hơn và đẹp hơn. Các input về dữ liệu và các biên ràng buộc, được chương trình thiết kế rất thân thiện với người sử dụng. 2. HỌ SPLINE VỚI BIÊN TỰ NHIÊN Ở đây ta sẽ xây dựng các hàm spline cấp 3 (đa thức cấp 3) trong từng đoạn thoả mãn các yêu cầu trên qua định nghĩa sau: Định nghĩa: Hàm  f x xác định trên  ,a b , một phân hoạch của  ,a b là bộ n+1 điểm: 1 2 3 1... n na x x x x x b       . Một spline đa thức bậc ba g(x) nội suy f(x) phải thoả mãn các điều kiện sau: (*) TS, Khoa Cơ bản, Trường Đại học Công nghiệp TP.HCM a) g(x) có đạo hàm đến cấp 2 liên tục trên  ,a b b) Trên mỗi đoạn con  1,k kx x  , 1,k n  ,    kg x g x là một đa thức bậc 3. c)    k k kg x f x y  , 1, 1k n   . d) Thoả mãn một trong 2 điều kiện biên sau đây: i.    // //1 1 0ng x g x   (điều kiện biên tự nhiên) ii.    / /1 1g x f x ,     / / 1 1n ng x f x  (điều kiện biên ràng buộc) Người ta chứng minh được rằng nếu hàm f(x) xác định trên tập điểm nút   1 1 k n k k x    thì tồn tại một hàm spline duy nhất thoả bốn điều kiện trên (xem [8]). Cách xây dựng hàm spline xin xem ([7], [8]). Ở đây ta chỉ trưng ra thuật toán để xây dựng chương trình. Thuật toán spline với biên tự nhiên: Bước 1: Tính 1k k kh x x  , 1,k n  , và giải hệ phương trình với (n+1) ẩn 1 2 1, , ..., nm m m  , nhưng thực chất là (n-1) phương trình và (n-1) ẩn, vì điều kiện biên tự nhiên: 1 1 0nm m   1 1 1 1 1 1 1 1 1 0 6 3 6 n k k k k k k k k k k k k k m m h h h h y y y y m m m h h                    , 2,k n  , (1.1) Đây là phương trình Ax=b, có ma trận A là 3 đường chéo chính sau: 1 2 3 4 n 1n                  giải hệ này để tìm ra 2 3, , ..., nm m m Bước 2: Tính ,k kA B : 2 6 k k k k h A y m  , 2 1 1 6 k k k k h B y m   , với 1,k n  (1.2) Bước 3: Xây dựng spline  kg x trên mỗi đoạn 1,k kx x    0 1 2 3 4 5 6 7 8 9 10 -6 -4 -2 0 2 4 6 8       3 3 1 1 1 1 6 6 , , 1, k k k k k k k k k k k k k k k x x x x x x x x g x m m A B h h h h x x x k n                      (1.3) Ví dụ 1: Cho mốc nội suy: 0 2 3 5 8 10 2 3 4 5 3 6 x y  Ở đây có 6 mốc nội suy, vậy ta phải xây dựng 5 hàm spline trên 5 đoạn. Xây dựng ma trận 3 đường chéo ở vế trái (VT), vế phải (VP) của (1.1): VT của (1.1) VP của (1.1) 1.000000 0.166667 0.000000 0.000000 9.5000 0.166667 1.000000 0.333333 0.000000 -6.5000 0.000000 0.333333 1.666667 0.500000 -1.1667 0.000000 0.000000 0.500000 1.666667 2.1667 Giải hệ phương trình trên tìm km , Tính các hệ số kA , kB từ (1.2) ta có: km , 1,6k  kA , 1,5k  kB , 1,5k  0 2.0000 -10.2831 10.9246 -4.8208 5.4246 -8.5477 9.6985 4.5461 0.6808 3.9788 1.3564 1.0958 2.2695 6.0000 0 Chương trình xuất ra 5 biểu thức hình thức hàm spline trên từng đoạn, được vẽ như (đồ thị 1): g1(t) = .91038*t^3+2-6.1415*t g2(t) = -1.4246*(t-2)^3+1.8208*(3-t)^3-25.3116+10.2454*t g3(t) = .56735e-1*(t-3)^3-.71231*(5-t)^3+17.4267-2.5761*t g4(t) = .60875e-1*(t-5)^3+.37823e-1*(8-t)^3+8.34980-.87418*t g5(t) = .91313e-1*(10-t)^3-12.6530+1.8653*t Đồ thị 1: Biểu diễn 5 hàm g1(t), g2(t), g3(t), g4(t), g5(t) trên 5 đoạn tương ứng. Nếu viết dưới dạng triển khai phân số thì hàm gồm 5 nhánh sau :                           2 108521 17670 x 32173 35340 x3 x 2    20760 589 989711 17670 x 29373 1178 x2 28673 8835 x3 x 3     86163 1178 462767 8835 x 6594 589 x2 13589 17670 x3 x 5    106579 5301 63091 17670 x 1 186 x2 611 26505 x3 x 8    46331 589 451091 17670 x 3227 1178 x2 3227 35340 x3 otherwise Ví dụ 2: Cho mốc nội suy: 0 2 3 5 6 8 11 13 3 4 4 2.5 2.5 3.5 3 3.5 x y Dùng 7 hàm spline xấp xỉ, đồ thị 3 g1(t) = -.27372e-1*t^3+3.0+.6095*t g2(t) = -.17154*(t-2)^3-.54744e-1*(3-t)^3+3.8211+.1168*t g3(t) = .83491e-1*(t-3)^3-.85768e-1*(5-t)^3+8.96738-1.42706*t g4(t) = .91182e-1*(t-5)^3+.16698*(6-t)^3+1.9540+.758e-1*t (1.4) g5(t) = -.53518e-1*(t-6)^3+.45591e-1*(8-t)^3-3.2438+.8965*t g6(t) = .24593e-1*(t-8)^3-.35679e-1*(11-t)^3+10.13644-.70913*t g7(t) = .36889e-1*(13-t)^3-1.6688+.3976*t 1x 2x nx 1nx   1g x   ng x  x 0 2 4 6 8 10 12 14 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 Đồ thị 2: Biểu diễn 7 hàm g1(t), g2(t), g3(t), g4(t), g5(t) , g6(t), g7(t) trên 7 đoạn tương ứng. 3. HỌ SPLINE VỚI BIÊN RÀNG BUỘC Điều kiện biên ràng buộc ở hai đầu mốc nội suy là: (,  cho trước)    / /1 1 1g x f x   ,     / / 1 1n n ng x f x    .  k=1, Trong phần chứng minh [7] đã có:  / 2 k k k k k k k k h A B g x m h h     , ==> k=1,  / 1 1 11 1 1 1 12 h A B g x m h h     Và    1 1 1 6 k k k k k k k k k y y h B A m m h h        vậy    / 1 2 1 11 1 1 1 2 12 6 h y y h g x m m m h         , hay 1 1 2 11 2 13 6 h h y y m m h        k=n+1, Trong phần chứng minh [7] đã có:  / 1 1 11 1 12 k k k k k k k k h A B g x m h h          ==> k=n+1,  / 1 1 2 n n n n n n n n h A B g x m h h      vậy    / 11 1 1 2 6 n n n n n n n n n n h y y h g x m m m h          , hay 11 6 3 n n n n n n n h h y y m m h       . Thuật toán spline với biên ràng buộc: Bước 1: Tính 1k k kh x x  , 1, 1k n   , và giải hệ phương trình với (n+1) ẩn 1 2 1, , ..., nm m m  , 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1: 3 6 2, : 6 3 6 1: 6 3 k k k k k k k k k k k k k n n n n n n n h h y y k m m h h h h h y y y y k n m m m h h h h y y k n m m h                                       (3.1) Đây là phương trình Ax=b, có ma trận A là 3 đường chéo chính sau: 1 2 3 4 n 1n                      giải hệ này để tìm ra 1 2 1, , ..., ,n nm m m m  . Bước 2: Tính ,k kA B : 2 6 k k k k h A y m  , 2 1 1 6 k k k k h B y m   , với 1,k n  (3.2) Bước 3: Xây dựng  kg x trên mỗi đoạn 1,k kx x    :       3 3 1 1 1 1 6 6 , , 1, k k k k k k k k k k k k k k k x x x x x x x x g x m m A B h h h h x x x k n                      (3.3) Ví dụ 3: Cho mốc nội suy: 0 2 3 5 6 8 11 13 3 4 4 2.5 2.5 3.5 3 3.5 x y Đây là ví dụ 3 ở trên, mục đích của ta là để so sánh 2 spline với 2 điều kiện biên tự nhiên và điều kiện biên ràng buộc. Với  =0, =0 cho trước:  =  / 1 0f x  , =   // 1 0nf x   Ở đây ta dùng 7 hàm spline xấp xỉ cho 7 đoạn tương ứng: Và chương trình xuất ra 7 biểu thức hàm như sau: g1(t) = -.59056e-1*t^3+.92028e-1*(2-t)^3+2.2638+1.1043*t g2(t) = -.15944*(t-2)^3-.11811*(3-t)^3+4.0355+.413e-1*t g3(t) = .81184e-1*(t-3)^3-.79719e-1*(5-t)^3+8.81872-1.39364*t g4(t) = .94663e-1*(t-5)^3+.16237*(6-t)^3+1.9991+.677e-1*t (3.4) g5(t) = -.57586e-1*(t-6)^3+.47331e-1*(8-t)^3-3.3962+.9196*t g6(t) = .32859e-1*(t-8)^3-.38391e-1*(11-t)^3+11.00004-.80793*t g7(t) = -.55894e-1*(t-11)^3+.49288e-1*(13-t)^3-4.7732+.6708*t Nếu viết dưới dạng triển khai phân số thì hàm gồm 7 nhánh sau : 3 2 1 202 471 3 1337 853 g x x     , 2x  3 2 2 126 498 839 1984 3049 4681 637 935 g x x x      , 3x  Bien tu nhien  // 1 0nf x   Bien rang buoc  / 1 0f x  3 2 3 563 969 1491 227 3499 503 220 68 g x x x    , 5x  3 2 4 165 2783 197 7319 2437 1852 19 290 g x x x      , 6x  3 2 5 431 1801 3827 3261 4108 829 266 98 g x x x      , 8x  3 2 6 227 1963 6453 6546 3186 955 332 115 g x x x    , 11x  3 2 7 613 1970 10349 15478 5828 523 232 87 g x x x      , 13x  Đồ thị 3: Trộn 7 biểu diễn 7 hàm gi(t), i=1..7, trên 7 đoạn tương ứng với biên ràng buộc của ví dụ 3 và 7 biểu diễn của 7 hàm của ví dụ 2 với biên tự nhiên. So sánh 7 biểu thức hàm (1.4) và (3.4), ta thấy rằng: hai nhánh g1(x) và g7(x) của ví dụ 3 và ví dụ 4 khác nhau nhiều, còn các nhánh còn lại bên trong vẫn khác nhau chút ít, do sự khác biệt ở hai đầu mút biên và chúng truyền vào bên trong. Và ở đồ thị 4 thể hiện rất rõ điều kiện biên ở hai đầu mốc nội suy là:    / /1 1 1 0g x f x    ,     / / 7 8 8 0g x f x    . 4. SO SÁNH SPLINE VÀ XẤP XỈ THEO NEWTON Phần này ta lấy lại ví dụ 4 với mốc nội suy: 0 2 3 5 6 8 11 13 3 4 4 2.5 2.5 3.5 3 3.5 x y thì đa thức xấp xỉ theo Newton (hay Lagrange cũng cho cùng kết quả) là đa thức bậc 7 và có biểu thức phức tạp như sau:   7 6 5 4 3 2 1 7 8 54 211 846 1853 851 1960 3 130071 18713 3966 1745 828 181 639 t t t t t t t P t          Xem đồ thị 5, ta thấy đa thức này xấp xỉ ở hai đầu mút nội suy là không chính xác, và chính vì thế nó gây ra sai số rất lớn. Đồ thị 4: So sánh giữa 2 đường cong xấp xỉ spline và Newton . 5. KẾT LUẬN Nhận xét họ hàm spline đẹp hơn và tốt hơn (theo nghĩa xấp xỉ ít sai số) so với cách dựng của hai phương pháp nội suy Lagrange, nội suy Newton, nhưng bù vào đó là tính toán phức tạp hơn. Chính vì thế ngày nay các họ spline được thay thế cho các đường cong xấp xỉ theo Newton hay Lagrange. TÀI LIỆU THAM KHẢO 1. Nguyễn Phú Vinh, Giáo trình phương pháp tính, Đại học Công nghiệp, 2006. 2. Nguyễn Phú Vinh, Các chương trình tính toán bằng matlab, Đại học Công nghiệp, 2007. Newton spline 3. Nguyễn Phú Vinh và các tác giả khác. Các giáo trình toán cao cấp đại học: A1, A2, A3. ĐHCN. 2006. 4. Nguyễn Phú Vinh, Các sources phương pháp tính, toán cơ được viết bằng VB for excel và Matlab, Computer Graphics. 5. Vera.B.Anand. Đồ họa máy tính và mô hình hoá hình học (Computer Graphics and geometric modeling). 6. Đêmiđôvich. Computatinal Mathematics. Mir. 1990. 7. Lê Thái Thanh, Lê Ngọc Lăng, Nguyễn Quốc Lân, Phương pháp tính. 2003. 8. Burden R.L, Faires T.D. Numerical Analysis, Brooks/Cole Publishing Company, USA- 1997.