Thuật toán từng bước (Có độ dốc) tốt hơn cho việc giải hệ phương trình phi tuyến

TAÏP CHÍ ÑAÏI HOÏC SAØI GOØN Soá 6 - Thaùng 6/2011 91 THUẬT TOÁN TỪNG BƯỚC (CÓ ĐỘ DỐC) TỐT HƠN CHO VIỆC GIẢI HỆ PHƯƠNG TRÌNH PHI TUYẾN NGUYỄN PHÚ VINH (*) TÓM TẮT Trong bài báo này chúng tôi trình bày phương pháp có tên “từng bước tốt hơn”: biến việc giải hệ phương trình phi tuyến thành bài toán tìm cực tiểu của một hàm có dạng tổng bình phương. Phương pháp này dựa vào hướng gradient giảm dần giá trị của hàm để dần xấp xỉ đến giá trị cực tiểu của hàm và giá trị đó cũng chính là nghiệm địa phương của hệ phương trình ban đầu. Quá trình này đã được lập trình trên Matlab để thử nghiệm, so sánh nghiệm và tốc độ hội tụ của phương pháp này với với phương pháp Newton. ABSTRACT In this article, we study the steepest - descent algorithmwhich is a transformation of the solving nonlinear equations system into finding the minimum of a multivariable function belonging to the form of sum of squares. This method is on the decent of gradient of the function. This function value is approached by the minimum of the function which is the local solution of the above nonlinear equations system. This process is programmed by Matlab in order to test the solution of this problem. It is used to compare the solution and the rate of convergence of this method and Newton method. 1. ĐẶT VẤN ĐỀ (*) Việc giải xấp xỉ gần đúng một hệ phương trình phi tuyến là việc mà nhiều tác giả trong và ngoài nước đã từ lâu vẫn quan tâm. Nhưng chưa có một phương pháp nào tổng quát để giải mọi hệ phương trình. Hiện nay phương pháp Newton có thể xem là phương pháp tốt nhất, tốc độ hội tụ nhanh nhất để giải hệ phi tuyến, nhưng vẫn không phải là phương pháp vạn năng cho (*) TS, Khoa Cơ bản, Đại học Công nghiệp Tp. Hồ Chí Minh mọi hệ phương trình. Trong bài báo này chúng tôi trình bày phương pháp có tên “từng bước tốt hơn”: biến việc giải hệ phương trình phi tuyến thành bài toán tìm cực tiểu của một hàm có dạng tổng bình phương. Phương pháp này dựa vào hướng gradient giảm dần của hàm để dần xấp xỉ đến cực tiểu của hàm và đó chính là nghiệm địa phương của hệ phương trình ban đầu. 2. PHƯƠNG PHÁP VÀ KẾT QUẢ NGHIÊN CỨU 2.1. Phương pháp luận về thuật toán töøng böôùc toát hôn (steepest) Đặt bài toán: hãy giải hệ phi tuyến đa biến sau: 92       1 1 2 n 2 1 2 n n 1 2 n f x , x ,..., x = 0 f x , x ,..., x = 0 F = 0 - - - - - - - - - - - f x , x ,..., x = 0        , trong đó       1 1 2 n 2 1 2 n n 1 2 n f x , x ,..., x f x , x ,..., x F = ....................... f x , x ,..., x              ta biến bài toán thành bài toán tìm cực tiểu hàm có dạng tổng bình phương sau:        1 2 2 2 2 n i 1 1 2 1 2 n n 1 2 n2 n g x , x ,.., x = f = f x , x ,.., x + f x , x ,.., x +..+ f x , x , .., xn 2 i=1  vậy bài toán đưa về là bài toán tìm p địa phương:     nx R g p = min g x  =0. Để làm giảm từng bước hàm g(x) và dần đến 0 là cực tiểu của hàm, tiến hành với thuật toán như sau: bắt đầu với (0)x , sau đó tìm điểm dần tốt hơn (1)x sao cho:  (1) (0) (0)x = x -α g x với  >0, (chú ý g là hướng có khả năng làm giảm hàm g(x)), ta sẽ điều chỉnh  như thế nào đó để:    (1) (0)g x g x . Để xác định  , cần tìm cực tiểu của hàm một biến     (0) (0)h g x g x    . Nhưng việc khảo sát và lấy đạo hàm trực tiếp của  h  một cách bài bản thì thật là khó, thay vào đó làm theo phương pháp steepest như sau: Chọn ba số 1 2 3    và hi vọng chúng sẽ gần vị trí cực tiểu của hàm  h  , bằng cách ta xây dựng một tam thức bậc hai P(x) chính là đa thức nội suy Newton tiến của  h  tại ba điểm mốc nội suy 1 2,  và 3 . Sau đó ta tìm  1 3,   để  P  là cực tiểu trong đoạn  1 3,  và dùng  P  để xấp xỉ cực tiểu của  h  . Rồi dùng  làm giá trị xấp xỉ mới cho hàm g, tức là:  (1) (0) (0)x x g x   . Tổng quát ta có công thức truy hồi để tìm nghiệm xấp xỉ:  ( 1) ( ) ( )i i ix x g x    Tính công thức gradient hàm g: 93       n n i i i i,1 i,2 i,3 i,n i=1 i=1 g x = 2f f = 2f f ,f ,f ,..,f =  1 1,1 1,2 1,n2f f , f ,.., f +  2 2,1 2,2 2,n2f f , f ,.., f +    3 3,1 3,2 3,n n n,1 n,2 n,n2f f , f ,.., f +...+ 2f f , f ,.., f = = T 1,1 2,1 n,1 1,1 1,2 1,n1 1 1,2 2,2 n,2 2,1 2,2 2,n2 2 n n1,n 2,n n,n n,1 n,2 n,n f f M f f f M ff f f f M f f f M ff f 2 = 2 M MM M M M M M M M f ff f M f f f M f                                              =  T 1 2 n2J F x , x ,..., x Trong đó J là Jacobi của hàm  1 2, , ..., nF x x x . Sau đây là thuật toán steepest. 2.2. Thuật toán từng bước tốt hơn (steepest) Nhập giá trị ban đầu:  1 2, , ..., T nx x x x , số bước N cần tính, sai số chịu đựng tính toán TOL (tolerant) while k N do 1. Set  1 1 2, , ..., ng g x x x    1 2 1 2, , ..., 2 , , ..., T n nz g x x x J F x x x   0 2 z z 2. If 0z 0 then OUTPUT(“Zero gradient, may have a minimum”) , STOP 3. 0 z z z  , 1 0  , 3 1  ,  3 3g g x z  4. While  3 1g g do steps 5 and 6 5. 3 3 / 2  ,  3 3g g x z  6. If 3 / 2TOL  then OUTPUT(“Improvement, may have a minimum”) , STOP 7. 2 3 / 2  ,  2 2g g x z  8. 2 11 2 1 g g h      , 3 22 3 2 g g h      , 2 13 3 1 h h h      9. 10 2 3 0.5 h h          ,  0 0g g x z  10. Tìm  từ  0 3,  để    0 3min ,g g x z g g   11. x x z  12. If 1g g TOL  then OUTPUT(“successful”) , STOP 94 13. k=k+1 End 2.3. Các ví dụ và so sánh kết quả Ví dụ 1: giải hệ:     1 2 1 1 2 3 2 2 2 1 3 3 3 6 2cos 1 0 9 sin 1.06 0.9 0 60 3 10 3 0 x x f x x x f x x x f x e                       , lấy nghiệm ban đầu  (0) 1,1,1 Tx  . Với hàm 2 2 2 21 2 3 1 n i i g f f f f      , Tính Jacobi: J=           1 2 1 2 2 3 3 2 3 2 3 1 2 2 1 3 1 3 2 1 6 2sin 2sin 2cos 9 sin 1.06 sin 1.06 3 3 60 x x x x x x x x x x xx x x x x x e x e                          , và 1 2 3 f F f f              k=1, Thế (0)x để tính  (0)1g g x  8163.7523; z= 2 TJ F = -136.9376876467 24.4584536238 10759.2205558461          , tính chuẩn 0 2 z z =10760.1197537; và  (0) 0 -0.01272640925754 0.00227306518733 0.99991643234926 g x z z           . Vậy Với    (0) (0)1 1 10 g g x z g x       8163.7523 Với 3 1  F  (0) 3x z = 3.07635846249709 11.32373711890548 29.51313389084277          95 ==>  (0)3 3g g x z  1008.71607578 Vì g3 < g1 nên ta chấp nhận 3 1  và gán 2 0.5  , rồi tính F  (0) 2x z = 3.28250948349129 11.48734095680031 59.51632652171600          ==>  (0)2 2g g x z   3684.9269934065 Bây giờ dựa vào bảng mốc nội suy: 1 2 3 1 2 3g g g    , để nội suy ra hàm bậc hai theo  :          1 1 1 3 1 2 1 1 3 2P g h h g h h                    bằng phương pháp Newton tiến nhờ bảng dữ liệu sau: 1 1 2 2 3 3 0 8163.7523 0.5 3684.9269 1 1008.7160 g g g          ; ta có: 2 11 2 1 g g h      =-8957.6507317, 3 22 3 2 g g h      = -5352.42183523, 2 13 3 1 h h h      =3605.228896495 Thế 1 3,h h vào  P  , sau đó ta lấy đạo hàm  P  để tìm cực tiểu của nó, tức giải  / 0P   để tìm được nghiệm: 10 2 3 0.5 h h            =1.4923137322, ==> F  (0) 0x z = 3.34977423120232 11.14453481099873 -0.02878932588947          ==>  (0)0 0g g x z  =135.4224723 ta thấy 0g nhỏ hơn 1g và 3g . Vậy ta lấy (1) (0) 0 1.01899179529672 0.99660787360675 -0.49218902305463 x x z             và tính  (1)g x  135.4224723 k=2, tính lại z= 2 TJ F = 58.072382205939 203.772128736778 -2.042514338131          .............  (2) (1) 0 0.695645, -0.137993, -0.480816 T x x z   và  (2)g x 10.107836...... 96 Và thuật toán cứ thế lặp lại với công thức truy hồi:  ( 1) ( ) ( )i i ix x g x    cho đến khi đạt được nghiệm khá gần với nghiệm chính xác là: (0.498145, -0.199606, -0.528826) T ở đây ta lấy tiêu chuẩn dừng của thuật toán là: ( ) ( 1) 5 10 k k x x      , khoảng 28 bước (xem bảng sau) mới đạt được điều đó. k ( ) 1 k x ( ) 2 k x ( ) 3 k x g(x) 1 1.018992 0.996608 -0.492189 135.422472 2 0.695645 -0.137993 -0.480816 10.107836 3 0.693185 -0.137886 -0.528789 1.814149 4 0.583604 -0.226612 -0.523367 0.494123 5 0.582668 -0.225863 -0.530753 0.295721 6 0.555173 -0.209607 -0.525621 0.182806 7 0.554550 -0.209322 -0.529834 0.118411 8 0.535784 -0.204358 -0.526724 0.077100 9 0.535372 -0.204218 -0.529429 0.050569 1 0 0.522932 -0.202299 -0.527436 0.033205 1 1 0.522660 -0.202218 -0.529208 0.021829 1 2 0.514463 -0.201270 -0.527910 0.014356 1 3 0.514284 -0.201220 -0.529074 0.009448 1 4 0.508890 -0.200672 -0.528222 0.006219 1 5 0.508772 -0.200640 -0.528988 0.004094 1 6 0.505221 -0.200300 -0.528428 0.002696 1 7 0.505144 -0.200279 -0.528932 0.001775 1 8 0.502806 -0.200060 -0.528564 0.001169 1 9 0.502755 -0.200046 -0.528896 0.000770 2 0 0.501215 -0.199904 -0.528653 0.000507 97 2 1 0.501181 -0.199895 -0.528872 0.000334 2 2 0.500167 -0.199802 -0.528712 0.000220 2 3 0.500145 -0.199796 -0.528856 0.000145 2 4 0.499477 -0.199735 -0.528751 0.000096 2 5 0.499463 -0.199731 -0.528846 0.000063 2 6 0.499023 -0.199691 -0.528777 0.000041 2 7 0.499013 -0.199688 -0.528839 0.000027 2 8 0.498723 -0.199662 -0.528793 0.000018 Vậy sau 28 bước ta có nghiệm xấp xỉ là: (498145, -0.199606, -0.528826)T So với phương pháp Newton chính thống với ví dụ trên như bảng sau: ( ) 1 k x ( ) 2 k x ( ) 3 k x 1.127638 -0.270927 -0.513022 0.498513 -0.192263 -0.523877 0.498150 -0.199606 -0.528826 0.498145 -0.199606 -0.528826 Ta thấy phương pháp Newton chính thống có tốc độ hội tụ nhanh hơn nhiều so với phương pháp steepest. Ta xét thêm ví dụ sau: Ví dụ 2: Giải hệ  1 1 1 2 3 24 2 1 2 3 3 2 2 3 1 2 2 3 cos 1 0 1 0.05 0.15 1 0 0.1 0.01 1 0 f x x x x f x x x x f x x x x                        , lấy  (0) 0,0,0 Tx  ; Kết quả được chương trình in ra như bảng sau: chỉ cần năm bước là hội tụ k ( ) 1 k x ( ) 2 k x ( ) 3 k x g(x) 1 0.000000 0.009944 0.994385 0.008090 2 -0.020018 0.090056 0.995265 0.000449 3 -0.000128 0.094959 1.000282 0.000025 98 4 -0.001133 0.099438 0.999764 0.000001 5 -0.000003 0.099718 0.999995 0.000000 Nghiệm đúng của hệ này là: x1=0, x2=0.1, x3=1. 3. KẾT LUẬN Phương pháp steepest có tốc độ hội tụ chậm hơn nhiều so với phương pháp chính thống Newton, nhưng phương pháp steepest luôn chọn được hướng (gradient tốt dần) hội tụ của bài toán, thậm chí có những bài toán mà phương pháp Newton không cho ta kết quả, trong khi đó phương pháp steepest lại cho ta kết quả, tất nhiên tốc độ hội tụ của phương pháp này rất chậm. Trên đây cũng chỉ là một trong nhiều phương pháp mà thế giới đã nghiên cứu. Ở đây chúng tôi dựa vào phương pháp tối ưu để tìm cực trị của hàm phi tuyến đa biến cho trường hợp riêng khi hàm phi tuyến đa biến lại là dạng tổng bình phương của các hàm phi tuyến đa biến. Các kết quả này chỉ là bước đầu của việc tìm nghiệm địa phương, chưa đề cập đến việc đánh giá sai số giữa nghiệm xấp xỉ và nghiệm chính xác, cũng như chưa xác định vùng cận để chọn giá trị ban đầu (0)x sao cho bài toán hội tụ đến nghiệm địa phương cần tìm, chúng tôi sẽ hẹn những bài báo sau.