Phương pháp Newton-Raphson

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 \mathbf{Ax}=\mathbf{b} (với \mathbf{A} là ma trận n \times 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 f(x), mục tiêu là tìm giá trị của x thỏa phương trình f(x) = 0. (Tất nhiên trong trường hợp f(x) 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:

\displaystyle f'(x_1) = \frac{df(x_1)}{dx} \approx \frac{f(x_1) - f(x_2)}{x_1 - x_2}

\displaystyle \Rightarrow x_2 = x_1 - \frac{f(x_1)}{f'(x_1)} + \frac{f(x_2)}{f'(x_1)}

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

\displaystyle x_2 = x_1 - \frac{f(x_1)}{f'(x_1)}

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à:

\displaystyle x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

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

Hình 1: Phương pháp Newton-Raphson

Hình 1: Phương pháp Newton-Raphson

 

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 n càng lớn, giá trị của x_n càng gần về vị trí giao với trục hoành của đồ thị hàm f(x).

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 x_0. Giá trị x_0 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 f(x) = 0 vô nghiệm, hoặc khi giá trị đoán mò x_0 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 f(x) là đa thức theo x thì có thể dễ dàng cài đặt phương pháp này trong MATLAB như sau

function [ x ] = nrsolve( ce, res, maxLoop )
  [ceX, ceY] = size(ce);

  if (ceX > 1 && ceY == 1)
    ce = ce';
    [ceX, ceY] = size(ce);
  elseif (~(ceX == 1 && ceY > 1))
    disp(' ce is not valid ')
    return;
  end
  res = abs(res);

  % calculate derivate coefficients
  dce = zeros(ceX, ceY  - 1);
  for i=2:ceY
    if (ce(i) ~= 0)
      dce(i-1) = ce(i)*(i-1);
    end
  end

  % a rude guess, I guess :-)
  x0 = 1;
  x1 = x0 - evalf(ce, x0)/evalf(dce, x0);
  x0 = x1;
  iCnt = 0;
  while (abs(evalf(ce, x1)) > res && iCnt < maxLoop)
    x1 = x0 - evalf(ce, x0)/evalf(dce, x0);
    x0 = x1;
    iCnt = iCnt + 1;
  end

  x = x1;
end

function val = evalf(cs, x)
  val = cs(1);
  for i=2:size(cs, 2)
    if (cs(i) ~= 0)
      val = val + cs(i)*x^(i-1);
    end
  end
end

Trong chương trình trên, ce là vector chứa hệ số của đa thức. Chẳng hạn với đa thức f(x) = x^5 + 10x^3 - \frac{1}{2}x - 4 thì ce = \left[-4, \frac{-1}{2}, 0, 10, 1\right]. Tham số res 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 maxLoop 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 x^5 - 5 = 0 (cũng chính là tìm xấp xỉ giá trị \sqrt[5]{5}) như sau:

nrsolve([-5 0 0 0 0 1], 0.00000000001, 100)
ans =
    1.3797

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

>> tic
nrsolve([-1000 0 0 0 0 1], 0.00000000001, 100)
toc
tic
1000^(1/5)
toc

ans =
 3.9811
Elapsed time is 0.004110 seconds.

ans =
 3.9811
Elapsed time is 0.002021 seconds.

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

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

\displaystyle \begin{cases} f_1\left(x_1, x_2, ..., x_n\right) = 0 \\ f_2\left(x_1, x_2, ..., x_n\right) = 0 \\ ... \\ f_n\left(x_1, x_2, ..., x_n\right) = 0 \end{cases}

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 f_j:

\displaystyle df_j = \frac{\partial f_j}{\partial x_1}dx_1 + ... + \frac{\partial f_j}{\partial x_n}dx_n = \sum_{i=1}^n\frac{\partial f_j}{\partial x_1}dx_1

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

\displaystyle f_j\left(\mathbf{x}^{(2)}\right) - f_j\left(\mathbf{x}^{(1)}\right) = \sum_{i=1}^n\left(\frac{\partial f_j}{\partial x_i}\right)\left(x_i^{(2)}-x_i^{(1)}\right)              (1)

với \mathbf{x}^{(k)} là giá trị của nghiệm ở lần lặp thứ k. Ta mong muốn \displaystyle f_j\left(\mathbf{x}^{(2)}\right) = 0, do đó hệ n phương trình (1) có thể viết dưới dạng ma trận:

\displaystyle J^{(k)}\mathbf{\delta x}^{(k)} = -R^{(k)}      (2)

trong đó J^{(k)} là ma trận Jacobi kích thước n\times n tại bước lặp thứ k

\displaystyle \left(J^{(k)}\right)_{ji} = \left(\frac{\partial f_j}{\partial x_i}\right)      (3)

R^{(k)} là vector thặng dư kích thước n \times 1 tại bước lặp thứ k:

\displaystyle \left(R^{(k)}\right)_j=f_j\left(\mathbf{x}^{(k)}\right)        (4)

và cuối cùng

\displaystyle \mathbf{\delta x}^{(k)} = \mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}            (5)

Do đó nghiệm của hệ tại bước thứ k+1

\displaystyle \mathbf{x}^{(k+1)} = \mathbf{\delta x}^{(k)} +\mathbf{x}^{(k)}         (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 \mathbf{x_o}

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ính

4. Cập nhật giá trị \mathbf{x} 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 f(x_2) - f(x_1) với x_2 “đủ gần” với x_1). 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 \mathbf{Ax}=\mathbf{b} thì trong MATLAB được cài bằng phép chia trái \mathbf{x} = \mathbf{A\backslash b}, 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 \mathbf{Ax}=\mathbf{b} có nhiều phương pháp giải, trong đó có Conjugate gradient, Cholesky decomposition… Những phương pháp này sẽ là chủ đề của bài tiếp theo.

About these ads

10 Responses to “Phương pháp Newton-Raphson”

  1. Thiệt là thiệt thòi cho dân tin tụi mình. Giờ phải vỡ lòng lại phương pháp tính. Bọn điện tử bên BK học rất chắc phần này. Càng nghĩ càng tức chương trình, nhưng mà nếu học đủ như bọn nó thì liệu mình có lập trình koool được ko :)))

  2. Tức cũng vậy thôi a :)) Cứ AQ là “mình chưa biết là do mình chưa được học”, vậy đi :D

  3. nếu đây là phương trình phi tuyến thì nhập vô làm sao?

  4. Hi bạn,
    chương trình MATLAB viết trong bài chỉ là để minh họa ý tưởng, do đó giả sử hàm f(x) là đa thức. Trong trường hợp tổng quát, do việc tính đạo hàm gặp khó khăn nên việc cài đặt sẽ phức tạp hơn, bạn có thể sử dụng các hàm fzero(), fsolve() của matlab.

  5. bạn nào giả được phuong trình phi tuýen này không, 485.x^0,025 + 200.x + 323,67, liên hệ với mình qua địa chỉ hoangmda80@gmail.com. cảm ơn các bạn nhiều.

  6. điện tử viễn thông cũng không chắc phương pháp tính bằng dân toán tin đâu,còn lập trình thì nói thật là CNTT chỉ là “thùng rỗng kêu to”,nếu nói là trong đầu toán ko bằng điện tử thì đừng nói đến lập tình.Bill Gate không bao giờ tuyễn nhân viên lập trình mà không giỏi toán đâu.:D.

  7. Anh ơi, bài viết của anh rõ ràng dễ hiểu ghê, nhưng mà em là dân Quant Finance nên cũng có nhìu cái hơi lơ mơ, cho phép em post question sau này nhé. Cảm ơn anh ^_^

  8. em đây học xây dựng mà họ cũng bắt hoc môn này khó quá anh nhem ak

Trackbacks

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 61 other followers

%d bloggers like this: