Phương pháp Newton-Raphson

Lê Văn Tuấn

Phương pháp Newton-Raphson

(Nguồn: https://phvu.net)

Trong chương trình phổ thông, khoảng lớp 5 học sinh được học cách giải phương trình bậc nhất 1 ẩn (qua một lớp bài toán được biết với cái tên nổi tiếng là tìm x). Khoảng cấp 2, HS được học cách giải phương trình bậc hai một ẩn và hệ phương trình bậc nhất hai ẩn.

Tốt nghiệp cấp 3, nếu đi học tiếp Cử nhân Công nghệ thông tin (như người viết bài này), HS sẽ được học tiếp cách giải hệ phương trình tuyến tính nhiều ẩn dạng  (với  là ma trận ) trong môn Đại số tuyến tính. Giải thuật lúc đó được dạy thì người viết không còn nhớ rõ, nhưng đại khái là tìm cách sử dụng các phép biến đổi dòng trên ma trận một cách khéo léo để tìm nghiệm. Mặc dù phương pháp đó có thể cài đặt được thành thuật toán, nhưng vì dựa vào trick (và “quan sát” theo kiểu của con người) khá nhiều nên không thực sự straightforward (sau vài ngày thì người viết nhớ ra tên phương pháp là Phép khử Gaussian)

Người viết có dịp tìm hiểu phương pháp Newton để giải hệ phương trình phi tuyến lần đầu tiên trong quá trình thực hiện luận văn. Sau đó gần đây lại có cơ hội làm việc với các phương pháp khác như Conjugate Gradient, LSQR …

Một điều đáng ngạc nhiên là, theo quá trình học của bản thân người viết, những phương pháp giải xấp xỉ/chính xác các phương trình/hệ phương trình tuyến tính/phi tuyến đã không được hệ thống hoá một cách đầy đủ, cả trong chương trình Toán phổ thông lẫn Đại học không chuyên ngành Toán (mặc dù tư tưởng chung của các phương pháp này khá tự nhiên). Vì vậy chuỗi bài này muốn hệ thống lại các phương pháp trên. Bài đầu tiên sẽ nói về phương pháp Newton-Raphson. Các thuật toán khác và một số vấn đề liên quan sẽ được trình bày trong các bài sau. Các phương pháp này là bước đi ban đầu vào optimization, lĩnh vực mà người viết mong muốn được tìm hiểu sâu hơn.

Phương pháp Newton-Raphson giải một (hệ) phương trình thực bằng cách xấp xỉ dần giá trị của nghiệm thông qua một dãy số (hội tụ về giá trị thật của nghiệm?). Chúng ta sẽ bắt đầu với trường hợp hàm một biến, sau đó mở rộng cho trường hợp nhiều biến.

1. Trường hợp một biến

Cho hàm một biến , mục tiêu là tìm giá trị của  thỏa phương trình . (Tất nhiên trong trường hợp  là hàm bậc nhất theo x thì không còn gì phải bàn, ở đây ta muốn xét đến trường hợp tổng quát).

Phương pháp Newton xuất phát từ việc xấp xỉ giá trị của đạo hàm tại một điểm:

Do ta đang muốn tìm nghiệm  sao cho , do đó với  thì phương trình trên trở thành:

Tổng quát, phương pháp Newton lặp liên tục để tìm nghiệm ngày càng gần giá trị thực, với công thức nghiệm là:

với giá trị khởi tạo  nào đó.

Hình 1 cho một cái nhìn trực quan hơn về lời giải của phương pháp này. Với  càng lớn, giá trị của  càng gần về vị trí giao với trục hoành của đồ thị hàm .

Tuy nhiên cũng giống như nhiều phương pháp tính xấp xỉ khác, khả năng tốc độ hội tụ của phương pháp Newton phụ thuộc nhiều vào giá trị “đoán mò” ban đầu . Giá trị  càng gần nghiệm thực của phương trình thì phương pháp này càng hội tụ nhanh chóng. Trong trường hợp xấu (khi  vô nghiệm, hoặc khi giá trị đoán mò  không cùng nằm trên một đoạn biến thiên của hàm số) thì phương pháp sẽ khó hội tụ nhanh. Trong trường hợp phương trình không có nghiệm, phương pháp này cũng không thể phát hiện.

Một điểm hạn chế nữa là do yêu cầu tính tích phân của hàm số, phương pháp này sẽ khó khăn trong trường hợp việc tính đạo hàm trở nên phức tạp. Trong trường hợp này, người ta có thể xấp xỉ việc tính đạo hàm bằng cách tìm phương trình đường thẳng đi qua hai điểm lân cận của hàm số tại vị trí đang xét (để xấp xỉ hệ số góc của tiếp tuyến – vốn có thể tính chính xác bằng đạo hàm tại điểm đó).

Do sự phức tạp trong việc tính đạo hàm nên việc cài đặt một chương trình tổng quát cho thuật toán Newton-Raphson là tương đối khó. Tuy nhiên trong trường hợp là đa thức theo  thì có thể dễ dàng cài đặt phương pháp này trong MATLAB như sau

Trong chương trình trên,  là vector chứa hệ số của đa thức. Chẳng hạn với đa thức  thì . Tham số là độ lỗi cho phép, chương trình dừng khi sai số nhỏ hơn ngưỡng này. Trong trường hợp không tìm được nghiệm, chương trình dừng sau  lần lặp.

Sử dụng hàm này, ta có thể tính nhanh nghiệm của phương trình  (cũng chính là tìm xấp xỉ giá trị ) như sau:

Một điều khá thú vị là sử dụng cặp hàm , ta thấy hàm  chạy không chậm hơn phép lũy thừa trong MATLAB nhiều lắm khi tính  (xét tới việc ta đã “đoán mò” giá trị  rất thô thiển).

2. Trường hợp nhiều biến

Giả sử ta phải giải hệ  phương trình  biến sau:

Một trong những phương pháp giải hệ phương trình này là phương pháp Newton nhiều biến (Multivariate Newton-Raphson Method – MNRM). Mở rộng từ phương pháp Newton cho hàm một biến, ta cũng bắt đầu từ việc đạo hàm từng phần mỗi phương trình :

Giống như trước, ta xấp xỉ như sau:

              (1)

với  là giá trị của nghiệm ở lần lặp thứ k. Ta mong muốn , do đó hệ  phương trình (1) có thể viết dưới dạng ma trận:

      (2)

trong đó  là ma trận Jacobi kích thước  tại bước lặp thứ 

      (3)

 là vector thặng dư kích thước  tại bước lặp thứ :

        (4)

và cuối cùng

            (5)

Do đó nghiệm của hệ tại bước thứ  là

         (6)

Sử dụng những phương trình này, thuật toán Newton tìm lời giải cho hệ phương trình rất đơn giản:

1. Đoán mò nghiệm 2. Tính ma trận Jacobi và vector thặng dư.3. Giải phương trình (2) bằng các phương pháp trong đại số tuyến tính4. Cập nhật giá trị  theo (6).5. Lặp lại bước 2 nếu lời giải chưa hội tụ.

3. Tổng kết

Phương pháp Newton-Raphson trong bài là một phương pháp hết sức đơn giản. Nó nằm trong họ các phương pháp Household mà nếu có dịp sẽ quay trở lại sau. Do hạn chế trong việc phải tính được đạo hàm của hàm số, phương pháp này thực sự có ý nghĩa khi biết trước dạng hàm. Trong trường hợp không biết trước dạng hàm thì cũng có thể áp dụng phương pháp này bằng cách tính xấp xỉ giá trị đạo hàm tại từng điểm (một trong những cách xấp xỉ là tính hiệu  với  “đủ gần” với ). Phương pháp xấp xỉ đạo hàm này đặc biệt có ý nghĩa khi cài đặt trên máy tính.

Trong MATLAB, việc tìm nghiệm của hệ các phương trình được cài đặt sẵn trong hàm fslove(), fzero()… Các hàm này cài đặt một số thuật toán khác để không phải tính đạo hàm như thuật toán Newton, hoặc tìm cách xấp xỉ giá trị đạo hàm. Đặc biệt hàm fsolve() không cần giá trị đoán mò ban đầu.

Một trường hợp đơn giản là giải hệ phương trình tuyến tính bậc nhất  thì trong MATLAB được cài bằng phép chia trái , một số thư viện số học khác như LAPACK, gsl, BLAS cũng có sẵn hàm cho trường hợp này. Trường hợp  có nhiều phương pháp giải, trong đó có Conjugate gradient, Cholesky decomposition